This R Markdown document is a final output for Tree Phenology Analysis using R (Winter Semester 2022). I am completing this file while sitting on my hard chair during a cold night in March. This document will be a “documentation” of how I learned R for the first time truly. With the help of the course modules, notes from the lectures, and of Mr. Google.
I hope this work will be of use somehow to the course’s creators– a humble perspective from an R beginner who strove to learn things from bottom upwards. For sure, I will be forever grateful to the existence of this course because it introduced me to the tool I hope I will be using often in my career. Prost!
Additionally, since I am a self-proclaimed nerd who likes thoroughly dissecting concepts, some chapters will start with condensed descriptions of the themes discussed in class (and I use mathematical and logic operation symbols often to show relationships), but these are not complete regurgitation of what was discussed; some are seemingly trivial terms which are new to me, followed by the important bit– the exercises and my take on them.
I hope this work of an R newbie is viewed with kindness and consideration. Admittedly, most of the time I just followed step by step the templates given as examples, and was not able to play around much with R features yet. As I came across problems, I attempted to troubleshoot and improvise what works and what does not, in order to progress in the exercises. My main aim in taking this course is to learn how R works in the context of scientific research, and in the end I think I achieved it somehow and more.
Dormancy prevents injury during winter. Triggered by environmental signals [photo-periods f(latitudes) + temperature f(geographical location, photoperiods, climate)], works through interaction of physiological processes.
Stages of dormancy:
Dormancy establishment- environment signals (declining temp and photoperiod)
Endo dormancy- tree’s endogenous properties dictate this (chilling temp)
Eco dormancy-acclimatization but not deeply dormant, no morphological /development changes. temperature- main environmental driver
(warm temperatures)
growth resumption
Physiological processes regulating dormancy:
Transport of water and solute stops (transport)
plant hormones and signaling compounds (Phyto-hormones)
actions of genes (genetics)
dynamics of non structural carbohydrates
Dormancy determination
Dormancy-> prevents injuries. Cold is required to regain the capabilities to grow
Endo-dormancy, eco-dormancy (growth is prevented in unfavorable condition, heat exposure)
Empirical-> growth chambers, shoots with buds collected in winter. 7 to 12 days. Buds growth or not. Longer exposure to chilling? Chill period, forcing period
Statistical-> Long phenological datasets ) temp records. Based of long series of dates, previous temperature records. Calculates max, mean temperatures, flowering dates. Chilling and heat accumulation. Flowering dates with mean and max temperatures. Number of independent variables> dependent variables (for PLS)
Chilling period-> when model coefficients are + significant values (temperature delays bloom dates)
Forcing period-> when model coefficients are -significant values (less temperature delays bloom dates)
experimental (empirical)
statistical- (dormancy determination) research purposes, collected in research centers. new commercial cultivars. Endodormancy break, chilling hours
Phenological datasets
Phenology in practical agriculture-> irrigation scheduling (when water is mostly taken up)
Tree-> budding, flowering, fruiting stages
Old vs BBCH scales
letters from A to J- every 2 days during spring. earliest and latest state
limitations; need for common phenology across the world
-need for standardized scale
-use of numbers is more practical than letters
Scales for plant growth stages
Cereals (BBCH, monocot, dicot)
Exercise:
Put yourself in the place of a breeder who wants to calculate the temperature requirements of a newly released cultivar. Which method will you use to calculate the chilling and forcing periods? Please justify your answer.
To calculate temperature requirements of a newly-released cultivar, I would do empirical first (trial and error, actual field experiments); but this could take time and could be expensive because the researcher would have to wait seasons by seasons and replicate to have enough repetitions, simulate the conditions through controlled environment (growth chamber). Statistical can be done with available historical data (long phenological data), and it has to be “applicable” to the new cultivar, meaning, more similarities than differences. However, the two methods could be done in parallel since they require different sets of data anyway. Results can be compared to verify if the models generated through statistical method is acceptable enough.
Which are the advantages (2) of the BBCH scale compared with earlier scales?
a. Standardized across the globe
b. More specific (has primary and secondary classifications in contrast to old version that used single letters)
Classify the following phenological stages of sweet cherry according to the BBCH scale:
budding- 54
flowering- 65
fruiting- 87
Drivers of climate change
Sun-cyclic variability (evidence: sunspots). Small portion
Aerosols- suspension of particles (dust, fires, seasalts, manmade), deflects some sunlight from Earth-> cooling effect. Major driver in industrial areas
Clouds- can have warming or cooling effect; depends on the type of cloud (varied effects depending on altitude, etc.). Lower= more sunlight reflected. Effect is very complex.
Ozone- tropospheric ozone (bad ozone, smog); stratospheric ozone (good ozone, blocks UV-B). +GHG -> warming effect. Destroyed by CFCs (ozone hole in S hemisphere, was fixed when the world took notice and acted-> Montreal protocol of 1987)
Surface albedo- properties of reflecting surface (land surface); light surfaces-> more reflection, dark surfaces-> less reflection
-deforestation- raises albedo
-drying-raises albedo
feedback loop-accelerates whatever is in place
-less ice/snow->lower albedo-> less heat reflection->warming
GHGs (CO2, CH4, N2O)- atmospheric gases that absorb long-wave radiation (traps heat)-> warming effect
-N2O, CH4~ agriculture
Long term drivers- trends in solar activity (star life cycle); ocean currents/continents; plants/animals (hypothetical)- imbalance of CO2 (heating due to higher CO2)-> flipping between “Snowball Earth” and “Greenhouse Earth”; volcanic and meteorite activity; Milankovic cycles (variation in tilt, eccentricity, precession of Earth’s orbit)
these factors are quantifiable: attribution studies
models have been crafted, compared to observed, and the recent increase(since industrialization) in global temperature can be explained by GHGE from human activities (burning of fossil fuels)
trends of CO2 emissions is still going up (new climatic territories)
degree of CO2 concentrations never seen by mankind, at high rates
Recent warming
Starting around 1980s, there is dramatic rise in global temperatures. Globally warmest years are very recent.
-> Siberia (2020) >8C warmer (possibility of +feedback=acceleration of warming effect)
Maunder minimum-> minimal sunspots, low solar radiation
Industrialization- GHGEs
Global surface temperatures started with warming, then a long period of cooling.
+1.5C is optimistic expectation (global temperature increase), but still requires adaptation-> +5C would cause major global changes unseen before by mankind
Rainfall is harder to predict than temperatures (no clear trends)
Moisture deficit in EU
Places with airports= good data (weather stations)
Future Scenarios
No way of validating which scenario is correct
GCM (Global Circulation Models)- computational effort (major drivers, feedbacks considered)
Climate drivers
Scenarios (fossil-intensive, low emissions, high emissions, intermediate)
RCPs (Representative Concentration Pathways)
2.6, 4.5, 6.0, 8.5- additional radiative forcing; higher number, more heating (higher emissions)
Downscaling (Dynamical downscaling, Statistical downscaling)- needed to increase resolution (more detail)
GCM>RCM>Impact Model
Statistical downscaling- higher resolution
-> use high resolution temp map, calculate diff between temperature map and average of pixel of RCM, delta change procedure (calculate bias), correct output using bias
in chillR, weather generator is implemented to reduce timestep of climate scenarios (temporal downscaling)
Climate change projection process
RCPs> GCM > downscaling (RCM) > statistical downscaling>Temperature, Precipitation, CO2 projections > Impact projection (e.g. what will happen to a certain species; biological response)
many ways to do process (no way to tell which is best), do all of the combinations and see what the models say
gap: in context of agriculture, models are created by data scientists who are not agricultural experts, there are gaps
Risk assessment (agricultural context)- additional step needed to project impacts in agriculture
Future temperatures
RCP 8.5: rise of 2 to 11C in global temperatures
Tipping point- point of sudden transition to new state (poorly considered in impact projections)
-> meltwater from Greenland dilutes salt concentration, less driver for global ocean currents (cold, salty water is heaviest); further slowing= transition of global climate to new state
-Greenland, arctic iceloss-> weakened albedo
-Permafrost thawing-> CO2 and methane release (more warming)
full of uncertainties
emissions can be translated to temperature increase
Paris Agreement: 1.5C increase only
-> 2060, CO2 neutral atleast
Impact Projection
system response
Statistical models- climate parameters x impact measure; used to explain past trends and project future impacts
-> species distribution modeling- climate parameters x presence/absence of species
ensemble modelling- using >1 approach, check overlaps
e.g. BiodiversityR package
limitations- things change, statistical relationships may not remain true; missing important factors; assumption: species distribution is in equilibrium with climate (x agricultural setups); must state uncertainties explicitly
Process-based models- computation, data, simulation to get quantitative projection of performance (e.g. crop models, phenology models, etc.); complex since need to include all relevant system components
projection of crop yields (limitations: does not consider real life threats-> diseases, weeds, etc.)
limitations- processes not completely understood; complex systems not modeled realistically; unclear uncertainties
complexity vs precision (more parameters, more errors)- aim for intermediate complexity to minimize errors “sweet spot”
express uncertainties and error estimates (always)
yield projections are complex, more uncertainties
before doing such projections, make sure you have a way to quantify the uncertainties!
Climate analogue models- comparing similar existing climates with projected climates; adaptation options
limitations: (analog vs target) other factors may be too similar-> cannot learn new things; other factors may be too different (e.g. cultural management, etc.)->too complex to model, huge uncertainties
Common limitations: climate data are scarce or poor quality; non-climatic factors also change through time, and difficult to project; CO2 can only be included in process-based models
Projecting climate-change impacts in complex agricultural systems:
need to understand what is happening: weather, climate, and crops + all other relevant factors affecting system performance; consider and communicate uncertainties
Exercise:
Main drivers of climate change (decade to century scale), explain mechanism through which the currently most important driver affects our climate.
at century scale: Sun (radiation), Aerosols (liquid,solid, mixed suspensions), Clouds (cooling and warming effect), Ozone (destroyed by CFCs, Montreal protocol of 1987), Surface albedo (reflects radiation) , GHGs, long term drivers (solar activity trends, ocean currents, plants vs animals, volcanic and meteorite activity, milankovic cycles)
mechanism through which most important driver affects climate: rising concentration of GHGs-> this is due to activities of mankind. CO2, CH4, and SO2 are released into the atmosphere. These stay in the atmosphere for very long periods and absorb long-wave radiation, which causes heating up of the Earth; thus causing “greenhouse effect”.
Explain briefly what is special about temperature dynamics of the recent decades, and why we have good reasons to be concerned.
Industrial revolution led to sharp rise in global temperatures (land, ocean surface temperatures) never seen before by mankind. A world that is quickly becoming >1.5C hotter than now would require tremendous adaption effort from future generations. Since we have not seen it before, the fear of uncertainty and looming feeling of doom is a great motivation to take GHGE reduction seriously; or else we might say ‘hi’ to dinosaurs sooner than we think.
What does the abbreviation ‘RCP’ stand for, how are RCPs defined, and what is their role in projecting future climates?
RCP: Representative Concentration Pathways. Each RCP is defined by additional radiative forcing; the higher the number, the higher the emissions and global heating. Currently, the commonly used ones are: RCP 2.6, 4.5, 6.0, 8.5, order is from the most conservative to the most emission-intensive estimate.
Impact projection starts at identifying RCPs -> inputs for climate change impact models. Process is as follows:
RCPs> GCM > downscaling (RCM) > statistical downscaling>Temperature, Precipitation, CO2 projections > Impact projection
Briefly describe the 4 climate impact projection methods described in the fourth video.
a. Statistical models- climate parameters x impact measure; used to explain past trends and project future impacts
b. Species distribution modelling- climate parameters x presence/absence of species.
c. Process-based models- use of mathematical equations to quantify system components and processes involved, and to project climate impact to the system. Complex since all relevant factors must be included (non-climatic); resource-intensive and requires processing of huge data
d. Climate analogue models- involves the use of a similar and existing location where the climate projection can be compared to. Target x analog locations
Chilling hours- hours of low temperature required to break fruit dormancy during winter. Various models have different definitions (Kaufmann, 2019):
Weinberger (1950): 1 CH=number of hours 0->7.2C
Utah model (1974): 1 CU= 1.4->12.4C
Dynamic model (1987): CP= -2-> 12C; optimal chilling 6->8C
Apples- 500 to 1000 CH
1 CP= 10 CH
Exercise:
Sketch out three data access and processing challenges that had to be overcome in order to produce chill projections with state-of-the-art methodology.
-limitation of equipment/available infrastructure in processing huge chunk of data
-incomplete or lack of high resolution climate datasets for some places; would need to make estimates (ground for errors) with respect to nearby places with available data
-form in which data is available in (especially for really large data); data in raster form adds additional steps before they can be rendered usable
Outline, in your understanding, the basic steps that are necessary to make such projections.
understand context and systems in location chosen (research + literature review)-> form your research question -> choose the appropriate chill model and climate projection ensembles to use -> obtain data needed for the chosen model-> interpolate/make estimates if there gaps -> process data using the model -> compare projections with existing projections/historic observations -> conclusion + expression of uncertainties and errors
CH computation requirement: hourly data;
read.table or read.csv functions;
.xls or .xlsx winters_hours_gaps dataset
kable format for improved aesthetic;
<- assigning new values to R objects comparison
commands (<, <=, ==,
=> and >);
c() command creates strings of numbers “vectors”;
& command combines comparisons;
FALSE and TRUE ~ 0 and
1 which command to identify items in dataset
of interest
Creating functions:
OurFunctionName <- function(argument1, argument2, ...) {ourCode}
Exercises:
Write a basic function that calculates warm hours (>25°C)
Year<-c(2000:2010)
Temps<-c(20:30)
randomdataset<-data.frame(Year, Temps)
##function creation, assigning of threshold to variable
WH<- function(randomdataset)
{
threshold_warm<- 25
randomdataset[,"Warm_year"]<- randomdataset$Temps>threshold_warm
return(randomdataset)
}
sampleWH<-WH(randomdataset)
write.csv(sampleWH, "data/sampleWH.csv", row.names = FALSE)
|
Year |
Temps |
Warm_year |
|---|---|---|
|
2000 |
20 |
FALSE |
|
2001 |
21 |
FALSE |
|
2002 |
22 |
FALSE |
|
2003 |
23 |
FALSE |
|
2004 |
24 |
FALSE |
|
2005 |
25 |
FALSE |
Apply this function to the Winters_hours_gaps dataset
library(knitr)
library(pander)
library(kableExtra)
library(chillR)
##function creation applied to Winter_hours_gaps; I made a new column that tells whether the temperatures are greater than the set threshold value: 25C
hourtemps<- Winters_hours_gaps[,c("Year", "Month", "Day", "Hour", "Temp")]
WH<- function(hourtemps)
{
threshold_warm<- 25
hourtemps[,"Warm_Hour"]<- hourtemps$Temp>threshold_warm
return(hourtemps)
}
write.csv(WH(hourtemps),"data/WH.csv", row.names = FALSE)
|
Year |
Month |
Day |
Hour |
Temp |
Warm_Hour |
|---|---|---|---|---|---|
|
2008 |
3 |
3 |
10 |
15.127 |
FALSE |
|
2008 |
3 |
3 |
11 |
17.153 |
FALSE |
|
2008 |
3 |
3 |
12 |
18.699 |
FALSE |
|
2008 |
3 |
3 |
13 |
18.699 |
FALSE |
|
2008 |
3 |
3 |
14 |
18.842 |
FALSE |
|
2008 |
3 |
3 |
15 |
19.508 |
FALSE |
Extend this function, so that it can take start and end dates as inputs and sums up warm hours between these dates
library(chillR)
##function creation applied to Winter_hours_gaps; I made a new column that tells whether the temperatures are greater than the set threshold value: 25C
hourtemps<- Winters_hours_gaps[,c("Year", "Month", "Day", "Hour", "Temp")]
WH<- function(hourtemps)
{
threshold_warm<- 25
hourtemps[,"Warm_Hour"]<- hourtemps$Temp>threshold_warm
return(hourtemps)
}
WH(hourtemps)[13:20,]
##sum of warm hours of start date and end year in YEARMODA
sum_WH<-function(hourtemps, Start_YEARMODA, End_YEARMODA)
{
Start_Year<-trunc(Start_YEARMODA/10000)
Start_Month<- trunc((Start_YEARMODA-Start_Year*10000)/100)
Start_Day<- Start_YEARMODA-Start_Year*10000-Start_Month*100
Start_Hour<-12
End_Year<-trunc(End_YEARMODA/10000)
End_Month<- trunc((End_YEARMODA-End_Year*10000)/100)
End_Day<- End_YEARMODA-End_Year*10000-End_Month*100
End_Hour<-12
Start_Date<-which(hourtemps$Year==Start_Year & hourtemps$Month==Start_Month & hourtemps$Day==Start_Day & hourtemps$Hour==Start_Hour)
End_Date<-which(hourtemps$Year==End_Year & hourtemps$Month==End_Month & hourtemps$Day==End_Day & hourtemps$Hour==End_Hour)
Warm_hours<- WH(hourtemps)
return(sum(Warm_hours$Warm_Hour[Start_Date:End_Date]))
}
##sample content for checking. # of warm hours in set time interval by manual counting= using sum_WH function), so I confirmed there's no error anywhere
sum_WH(hourtemps,20080517,20080518)
WH(hourtemps)[1803:1827,]
|
Year |
Month |
Day |
Hour |
Temp |
Warm_Hour |
|
|---|---|---|---|---|---|---|
|
1803 |
2008 |
5 |
17 |
12 |
31.663 |
TRUE |
|
1804 |
2008 |
5 |
17 |
13 |
32.639 |
TRUE |
|
1805 |
2008 |
5 |
17 |
14 |
33.183 |
TRUE |
|
1806 |
2008 |
5 |
17 |
15 |
33.417 |
TRUE |
|
1807 |
2008 |
5 |
17 |
16 |
31.791 |
TRUE |
|
1808 |
2008 |
5 |
17 |
17 |
29.790 |
TRUE |
|
1809 |
2008 |
5 |
17 |
18 |
29.790 |
TRUE |
|
1810 |
2008 |
5 |
17 |
19 |
27.087 |
TRUE |
|
1811 |
2008 |
5 |
17 |
20 |
28.196 |
TRUE |
|
1812 |
2008 |
5 |
17 |
21 |
28.493 |
TRUE |
|
1813 |
2008 |
5 |
17 |
22 |
27.875 |
TRUE |
|
1814 |
2008 |
5 |
17 |
23 |
21.843 |
FALSE |
|
1815 |
2008 |
5 |
18 |
0 |
19.175 |
FALSE |
|
1816 |
2008 |
5 |
18 |
1 |
18.176 |
FALSE |
|
1817 |
2008 |
5 |
18 |
2 |
18.533 |
FALSE |
|
1818 |
2008 |
5 |
18 |
3 |
16.201 |
FALSE |
|
1819 |
2008 |
5 |
18 |
4 |
15.008 |
FALSE |
|
1820 |
2008 |
5 |
18 |
5 |
14.721 |
FALSE |
|
1821 |
2008 |
5 |
18 |
6 |
15.199 |
FALSE |
|
1822 |
2008 |
5 |
18 |
7 |
17.938 |
FALSE |
|
1823 |
2008 |
5 |
18 |
8 |
23.088 |
FALSE |
|
1824 |
2008 |
5 |
18 |
9 |
25.162 |
TRUE |
|
1825 |
2008 |
5 |
18 |
10 |
27.161 |
TRUE |
|
1826 |
2008 |
5 |
18 |
11 |
29.615 |
TRUE |
|
1827 |
2008 |
5 |
18 |
12 |
31.459 |
TRUE |
##[1] 15CH model- not credible
Utah Model- more credible (different weights for different temperatures
Dynamic Model- complex, needs heavy maths.
Chilling_Hours() function determines if each hour are
within 0 to 7.2C ;
step_model() function- allows user to define model (set
thresholds and weights);
summ==TRUE cumulative sum;
data.frame- input, contains weights for temperature
ranges; custom() function implements chill model based on
set;
data.frame input make_JDay() function adds
Julian dates;
chilling() function implements CH, Utah, and Dynamic
models, and computes GDH (Growing Degree Days);
tempResponse function lets user pick what model they
want in chilling() function
Exercise:
Run the chilling() function on the
Winters_hours_gap dataset
library(chillR)
chill_sample<-chilling(make_JDay(Winters_hours_gaps),Start_JDay = 90, End_JDay = 100)
write.csv(chill_sample, "data/chill_sample.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Perc_complete |
Chilling_Hours |
Utah_Model |
Chill_portions |
GDH |
|---|---|---|---|---|---|---|---|---|
|
2007/2008 |
2008 |
11 |
11 |
100 |
40 |
15.5 |
2.009147 |
2406.52 |
Create your own temperature-weighting chill model using the
step_model() function
I looked at Milech et. al (2018) to see other existing models I have not seen yet so far in this module, and saw a nice table from their study summarizing models and the ranges. I picked Low Chill Model to be my guinea pig for this exercise portion; I did not use its equations, only the temperature ranges and weights.
##implementing the ranges from Low Chill Model
df<-data.frame(
lower= c(-1000, 1.8, 7.9, 13.9, 16.9, 19.4, 20.4),
upper=c(1.8, 7.9,13.9,16.9, 19.4, 20.4, 1000),
weight=c(0, 0.5, 1, 0.5, 0, -0.5, -1)
)
custom<- function(x) step_model(x, df)
##sample values for first 100 rows from the Winters_hours_gaps dataset)
custom(Winters_hours_gaps$Temp)[1:100]
write.csv(custom(Winters_hours_gaps$Temp), "data/custom.csv", row.names = FALSE)
Run this model on the Winters_hours_gaps dataset
using the tempResponse() function.
df<-data.frame(
lower= c(-1000, 1.8, 7.9, 13.9, 16.9, 19.4, 20.4),
upper=c(1.8, 7.9,13.9,16.9, 19.4, 20.4, 1000),
weight=c(0,0.5, 1, 0.5, 0, -0.5, -1)
)
custom<- function(x) step_model(x, df)
output<-tempResponse(make_JDay(Winters_hours_gaps), Start_JDay = 90, End_JDay = 100, models=list(custom=custom))
write.csv(output, "data/output.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Perc_complete |
custom |
|---|---|---|---|---|---|
|
2007/2008 |
2008 |
11 |
11 |
100 |
96.5 |
daylength function- generates daily temperature curves
from the latitude input;
ggplot2 package- for plotting;
JDay= Julian Date (1 to 365 or 366 in leap years);
melt command of the reshape2 package-
separates/melts dataframe, stacked; stack_hourly_temps()
function integrates daily dynamics using minimum and maximum temperature
data + latitude;
KA_weather dataset contains temperature data from Bonn
Uni; Empirical_daily_temperature_curve() function
determines typical temperatures at locations;
make_all_day_table function fills gaps in daily and hourly
temperature records
Exercise:
Choose a location of interest, find out its latitude and produce plots of daily sunrise, sunset and daylength.
I have always been fascinated with Nordic mythology– Vikings. Sometimes I ask myself: why don’t we just consult the three Norns who spin the webs of fate: Urd (Past), Verdandi (Present), and Skuld (future) instead of modelling climate scenarios? Would be a cooler way to predicting the heating of Earth! Below is a fun image of the three Norns I got from a Norse Mythology website.
My Nordic tale-loving soul convinced me to look at interesting locations in Norway, Sweden, and Denmark. I narrowed down my search by looking at the ranking of countries with the most number of weather stations, and voila! got my answer: Sweden. The plot below shows the number of weather stations across the globe from Jaffre (2019).
An interesting project implemented in Northern Sweden by Brännland Cider (a Swedish company making ice cider from apples) caught my attention. European innovation partnership for productivity and sustainability within agriculture (EIP-agri) of the Swedish Board of Agriculture granted them around 935,000 euro for this. The aim is to pilot climate-resilient agriculture in the Northern EU by trying to grow apples commercially in the harsh northern climate. The three counties involved are: Norrbotten, Västerbotten and Västernorrland; and I picked Västerbotten (64.34° N, 18.31° E) because it had many airports close to it– nearest was 12 km away. It would be very interesting to know what the field trials would look years in the future temperature-wise. Would they be able to make ice cider decades from now in that area?
I kept in mind that the location will be used repeatedly for the next topics so I chose quite strategically. Though my expectation is that they will be able to continue to grow apples because that region will not get warmer as per previous projections by Luedeling et al. (2011):
Now, for the exercise,
require(chillR)
require(ggplot2)
require(reshape2)
Days<-daylength(latitude=64.34, JDay=1:365)
Days_df<-data.frame(JDay=1:365, Sunrise=Days$Sunrise, Sunset=Days$Sunset, Daylength=Days$Daylength)
Days_df<-melt(Days_df, id=c("JDay"))
ggplot(Days_df, aes(JDay, value))+ geom_line(lwd=1.5) + facet_grid(cols=vars(variable))+ylab("Time of Day/Daylength (Hours)")+ theme_bw(base_size=20)
Produce an hourly dataset, based on idealized daily curves, for
the KA_weather dataset (included in
chillR)
require(chillR)
KA_weather_stacked<-stack_hourly_temps(KA_weather, latitude=64.34)
write.csv(KA_weather_stacked, "data/KA_weather_stacked.csv", row.names = FALSE)
|
hourtemps.Year |
hourtemps.Month |
hourtemps.Day |
hourtemps.Tmax |
hourtemps.Tmin |
hourtemps.JDay |
hourtemps.Hour |
hourtemps.Temp |
QC |
|---|---|---|---|---|---|---|---|---|
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
0 |
5.1 |
NA |
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
1 |
5.1 |
NA |
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
2 |
5.1 |
NA |
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
3 |
5.1 |
NA |
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
4 |
5.1 |
NA |
|
1998 |
1 |
1 |
8.2 |
5.1 |
1 |
5 |
5.1 |
NA |
Produce empirical temperature curve parameters for the
Winters_hours_gaps dataset, and use them to predict hourly
values from daily temperatures (this is very similar to the example
above, but please make sure you understand what’s going on)
library(chillR)
library(ggplot2)
##determining empirical daily temperatures from hourly data using Winters_hours_gaps dataset
coeffs<-Empirical_daily_temperature_curve(Winters_hours_gaps)
Winters_daily<-make_all_day_table(Winters_hours_gaps, input_timestep="hour")
Winters_hours<-Empirical_hourly_temperatures(Winters_daily,coeffs)
?Winters_hours_gaps
##input latitude of Winters, California (38.5 N) where Winters_hours_gaps dataset was recorded in
require(reshape2)
Winters_hours<-Winters_hours[,c("Year","Month","Day","Hour","Temp")]
colnames(Winters_hours)[ncol(Winters_hours)]<-"Temp_empirical"
Winters_ideal<-stack_hourly_temps(Winters_daily, latitude=38.5)$hourtemps
Winters_ideal<-Winters_ideal[,c("Year","Month","Day","Hour","Temp")]
colnames(Winters_ideal)[ncol(Winters_ideal)]<-"Temp_ideal"
##Triangular dataset
Winters_triangle<-Winters_daily
Winters_triangle[,"Hour"]<-0
Winters_triangle$Hour[nrow(Winters_triangle)]<-23
Winters_triangle[,"Temp"]<-0
Winters_triangle<-make_all_day_table(Winters_triangle,timestep="hour")
colnames(Winters_triangle)[ncol(Winters_triangle)]<-"Temp_triangular"
for(i in 2:nrow(Winters_triangle))
{if(is.na(Winters_triangle$Tmin[i])) Winters_triangle$Tmin[i]<-Winters_triangle$Tmin[i-1]
if(is.na(Winters_triangle$Tmax[i])) Winters_triangle$Tmax[i]<-Winters_triangle$Tmax[i-1]
}
Winters_triangle$Temp_triangular<-NA
Winters_triangle$Temp_triangular[which(Winters_triangle$Hour==6)]<-
Winters_triangle$Tmin[which(Winters_triangle$Hour==6)]
Winters_triangle$Temp_triangular[which(Winters_triangle$Hour==18)]<-
Winters_triangle$Tmax[which(Winters_triangle$Hour==18)]
Winters_triangle$Temp_triangular<-
interpolate_gaps(Winters_triangle$Temp_triangular)$interp
Winters_triangle<-Winters_triangle[,c("Year","Month","Day","Hour","Temp_triangular")]
##merging the dataframes
Winters_temps<-merge(Winters_hours_gaps,Winters_hours, by=c("Year","Month","Day","Hour"))
Winters_temps<-merge(Winters_temps,Winters_triangle, by=c("Year","Month","Day","Hour"))
Winters_temps<-merge(Winters_temps,Winters_ideal, by=c("Year","Month","Day","Hour"))
##conversion of dates, reorganizing dataframes, plotting
Winters_temps[,"DATE"]<-ISOdate(Winters_temps$Year,Winters_temps$Month, Winters_temps$Day, Winters_temps$Hour)
Winters_temps_to_plot<-Winters_temps[,c("DATE","Temp","Temp_empirical","Temp_triangular","Temp_ideal")]
Winters_temps_to_plot<-Winters_temps_to_plot[100:200,]
Winters_temps_to_plot<-melt(Winters_temps_to_plot, id=c("DATE"))
colnames(Winters_temps_to_plot)<-c("DATE","Method","Temperature")
ggplot(data=Winters_temps_to_plot, aes(DATE,Temperature, colour=Method)) +
geom_line(lwd=1.3) + ylab("Temperature (°C)") + xlab("Date")
GSOD- Global Summary of Day database
ChillR only has access to GSOD and California database
handle_gsod() function for data retrieval from GSOD;
Overlap_years column shows # years available;
Perc_interval_covered column shows % that target
(start-end date input) was covered
Exercise:
Choose a location of interest and find the 25 closest weather
stations using the handle_gsod function
Will use Västerbotten, Sweden (64.34° N, 18.31° E) again.
library(chillR)
##generation of station list based on target location's coordinates
station_list_Västerbotten<-handle_gsod(action="list_stations",
location=c(18.31,64.34),
time_interval=c(1990,2020))
write.csv(station_list_Västerbotten, "data/station_list_Västerbotten.csv")
|
chillR_code |
STATION.NAME |
CTRY |
Lat |
Long |
BEGIN |
END |
distance |
Overlap_years |
Perc_interval_covered |
|---|---|---|---|---|---|---|---|---|---|
|
022630_99999 |
FREDRIKA |
SW |
64.083 |
18.367 |
20010701 |
20230307 |
28.78 |
19.50 |
62.91 |
|
022640_99999 |
FREDRIKA |
SW |
64.083 |
18.417 |
19730101 |
20021028 |
29.12 |
12.83 |
41.37 |
|
023260_99999 |
FREDRIKA |
SW |
64.083 |
18.417 |
19770701 |
20041221 |
29.12 |
14.97 |
48.30 |
|
022455_99999 |
LYCKSELE |
SW |
64.548 |
18.716 |
19980408 |
20230307 |
30.33 |
22.73 |
73.33 |
|
022610_99999 |
LYCKSELE |
SW |
64.550 |
18.717 |
19730101 |
20230307 |
30.53 |
31.00 |
100.00 |
|
020630_99999 |
ASELE |
SW |
64.167 |
17.367 |
19550101 |
19770630 |
49.62 |
0.00 |
0.00 |
Download weather data for the most promising station on the list
From the 25 stations nearest the chosen location, the fifth (Lycksele, ChillR code 022610-99999) was the most promising; it covered 100% of the data for 1990 to 2020. Though it was quite far in terms of distance from the target area.
library(chillR)
##generation of station list based on target location's coordinates
station_list<-handle_gsod(action="list_stations",
location=c(18.31,64.34),
time_interval=c(1990,2020))
##having to enter just the number of station from the list created was a relief, thanks chillR coders!
weather<-handle_gsod(action="download_weather",
location=station_list$chillR_code[5],
time_interval=c(1990,2020))
##check sample data
weather[[1]]$weather[5090:6010,]
#downloading data
write.csv(weather,"data/Västerbotten_weather.csv",row.names=FALSE)
|
LYCKSELE.database |
LYCKSELE.weather.DATE |
LYCKSELE.weather.STN… |
LYCKSELE.weather.WBAN |
LYCKSELE.weather.YEAR |
LYCKSELE.weather.MONTH |
LYCKSELE.weather.DAY |
LYCKSELE.weather.TEMP |
LYCKSELE.weather.Count1 |
LYCKSELE.weather.DEWP |
LYCKSELE.weather.Count2 |
LYCKSELE.weather.SLP |
LYCKSELE.weather.Count3 |
LYCKSELE.weather.STP |
LYCKSELE.weather.Count4 |
LYCKSELE.weather.VISIB |
LYCKSELE.weather.Count5 |
LYCKSELE.weather.WDSP |
LYCKSELE.weather.Count6 |
LYCKSELE.weather.MXSPD |
LYCKSELE.weather.GUST |
LYCKSELE.weather.MAX |
LYCKSELE.weather.MaxFlag |
LYCKSELE.weather.MIN |
LYCKSELE.weather.MinFlag |
LYCKSELE.weather.PRCP |
LYCKSELE.weather.PrcpFlag |
LYCKSELE.weather.SNDP |
LYCKSELE.weather.FRSHTT |
LYCKSELE.weather.Year |
LYCKSELE.weather.Month |
LYCKSELE.weather.Day |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
GSOD |
1990-01-01 12:00:00 |
22610 |
99999 |
1990 |
1 |
1 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
1 |
|
GSOD |
1990-01-02 12:00:00 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
2 |
|
GSOD |
1990-01-03 12:00:00 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
3 |
|
GSOD |
1990-01-04 12:00:00 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
4 |
|
GSOD |
1990-01-05 12:00:00 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
5 |
|
GSOD |
1990-01-06 12:00:00 |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
NA |
1990 |
1 |
6 |
Convert the weather data into chillR format
library(chillR)
##generation of station list based on target location's coordinates
station_list<-handle_gsod(action="list_stations",
location=c(18.31,64.34),time_interval = c(1990,2020))
##having to enter just the number of station from the list created was a relief, thanks chillR coders!
weather<-handle_gsod(action="download_weather",
location=station_list$chillR_code[5],
time_interval=c(1990,2020))
##datacleaning, changing format to chillR
cleaned_weather<-handle_gsod(weather)
##check sample cleaned data, I used a non-boring range
cleaned_weather[[1]]$weather[5090:6010,]
##downloading data for future use
write.csv(cleaned_weather$weather,"data/Västerbotten_chillR_weather.csv",row.names=FALSE)
From the looks of the data I got as shown above from Västerbotten taken from GSOD database, it’s not possible to work with it for the next exercises. I already tried doing first exercise in Chapter 10 and it failed to be read using read.csv function (contains no data, that also happened to Bonn weather example from GSOD I ran). So for the next exercises, I would have to choose another location, possibly in Germany since the DWD database seems working well with chillR. To be able to move forward, I will also use DWD database for the filling of gaps exercise.
Another reason is we already know that Nordic regions will remain cold enough anyway as per previous projections. Makes me sad though; but will not cry because I am a tough Viking by heart.
linear interpolation can be done to fill in short gaps (hours to weeks), but longer gaps need better methods to reflect temperature dynamics
interpolate_gaps() function carries out simple
interpolation; fix_weather() function uses linear
interpolation between Tmax and Tmin columns and can be used to identify
gaps in the dataset
Exercise:
Use chillR functions to find out how many gaps you
have in this dataset (even if you have none, please still follow all
further steps)
I will repeat the previous steps done in the last exercise in this portion. New location of choice is the Swabian Orchard Paradise (48.4931° N, 9.3998° E) near Stuttgart. I had the chance to visit a Schwabian place before and ate some of their local dishes, courtesy of the tour by my mentor Martin Gummert. The orchard grows cherries, pears, and plums– and it is said that it has been centuries-old; would be interesting to check if they will survive another century.
A quick view of the place I got from Schwäbisches Streuobstparadies Website (n.d.):
Now onto the code:
library(chillR)
##generation of list, handle DWD dataset
station_list_dwd<-handle_dwd(action="list_stations",
location=c(9.40,48.49),
time_interval=c(19900101,20201231))
##check weather station list
require(kableExtra)
kable(station_list_dwd) %>%
kable_styling("striped", position = "left", font_size = 8)
##the list showed M<fc>nsingen-Apfelstetten (station 654) as the nearest and most number of years available (1946 to 2021), so I am choosing it (#2 in the station list)
##downloading the data from chosen weather station and checking data
weather_dwd<-handle_dwd(action="download_weather",
location=station_list_dwd$Station_ID[2],
time_interval=c(19900101,20201231))
head(weather_dwd$`M<fc>nsingen-Apfelstetten`$weather)
##datacleaning and checking data
cleaned_weather_dwd<-handle_dwd(weather_dwd)
cleaned_weather_dwd$`M<fc>nsingen-Apfelstetten`
##saving data in computer
write.csv(station_list_dwd,"data/station_list_dwd.csv",row.names=FALSE)
write.csv(weather_dwd,"data/MApfelstetten_weather1.csv",row.names=FALSE)
write.csv(cleaned_weather_dwd$`M<fc>nsingen-Apfelstetten`,"data/MApfelstetten_chillR_weather.csv",row.names=FALSE)
##assigning previous data to an object
MApfelstetten<-read.csv("data/MApfelstetten_chillR_weather.csv")
#checking for gaps using fix_weather function
MApfelstetten_QC<-fix_weather(MApfelstetten)$QC
write.csv(MApfelstetten_QC, "data/MApfelstetten_QC.csv", row.names = FALSE)
##saving and downloading the needed columns for next exercises
MApfelstetten_weather<-MApfelstetten[,c("Year","Month", "Day", "Tmax", "Tmin")]
kable(MApfelstetten_weather,) %>%
kable_styling("striped", position = "left", font_size = 10)
write.csv(MApfelstetten_weather, "data/MApfelstetten_weather.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Missing_Tmin |
Missing_Tmax |
Incomplete_days |
Perc_complete |
|---|---|---|---|---|---|---|---|
|
1989/1990 |
1990 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1990/1991 |
1991 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1991/1992 |
1992 |
366 |
366 |
0 |
0 |
0 |
100 |
|
1992/1993 |
1993 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1993/1994 |
1994 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1994/1995 |
1995 |
365 |
365 |
0 |
0 |
0 |
100 |
|
Year |
Month |
Day |
Tmax |
Tmin |
|---|---|---|---|---|
|
1990 |
1 |
1 |
0.6 |
-10.0 |
|
1990 |
1 |
2 |
4.3 |
-8.8 |
|
1990 |
1 |
3 |
-2.3 |
-4.7 |
|
1990 |
1 |
4 |
-3.3 |
-5.5 |
|
1990 |
1 |
5 |
0.4 |
-8.5 |
|
1990 |
1 |
6 |
0.5 |
-2.5 |
Surprise, surprise! according to the results from running the
fix_weather() function, there is no gap for the
Münsingen-Apfelstetten weather station (DWD database) for our
years of interest.
So now, I will create artificial gaps by deleting some data from the original dataset and work with filling it for the next exercises. The dataset I am using is here.
library(chillR)
##assigning previous data to an object
MApfelstetten_gaps<-read.csv("data/MApfelstetten_chillR_weather_gaps.csv")
#checking for gaps using fix_weather function
MApfelstetten_QC_gaps<-fix_weather(MApfelstetten_gaps)$QC
write.csv(MApfelstetten_QC_gaps, "data/MApfelstetten_QC_gaps.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Missing_Tmin |
Missing_Tmax |
Incomplete_days |
Perc_complete |
|---|---|---|---|---|---|---|---|
|
1989/1990 |
1990 |
365 |
365 |
38 |
38 |
38 |
89.6 |
|
1990/1991 |
1991 |
365 |
365 |
42 |
42 |
42 |
88.5 |
|
1991/1992 |
1992 |
366 |
366 |
0 |
0 |
0 |
100.0 |
|
1992/1993 |
1993 |
365 |
365 |
0 |
0 |
0 |
100.0 |
|
1993/1994 |
1994 |
365 |
365 |
0 |
0 |
0 |
100.0 |
|
1994/1995 |
1995 |
365 |
365 |
0 |
0 |
0 |
100.0 |
This confirmed the existence of gaps in 1990, 1991, 2004, 2005, 2012, 2018, 2019, and 2020.
Create a list of the 25 closest weather stations using the
handle_gsod function
I will use the location of the Swabian Orchard Paradise (48.4931° N,
9.3998° E) to check weather stations around via the
handle_dwd function.
library(chillR)
library(kableExtra)
##generation of list, handle DWD dataset
station_list_dwd<-handle_dwd(action="list_stations",
location=c(9.40,48.49),
time_interval=c(19900101,20201231))
|
Station_name |
Region |
Station_ID |
Elevation |
Latitude |
Longitude |
Begin |
End |
Distance |
|---|---|---|---|---|---|---|---|---|
|
Lenningen-Schopfloch |
Baden-W<fc>rttemberg |
2949 |
758 |
48.5396 |
9.5275 |
19470801 |
20020131 |
10.92 |
|
M<fc>nsingen-Apfelstetten |
Baden-W<fc>rttemberg |
3402 |
755 |
48.3851 |
9.4837 |
19460101 |
20211231 |
13.21 |
|
Reutlingen-Betzingen |
Baden-W<fc>rttemberg |
4166 |
360 |
48.5051 |
9.1833 |
19880501 |
20021231 |
16.10 |
|
N<fc>rtingen |
Baden-W<fc>rttemberg |
3671 |
280 |
48.6321 |
9.3261 |
19550101 |
19910331 |
16.72 |
|
Trochtelfingen-Wilsingen |
Baden-W<fc>rttemberg |
5107 |
769 |
48.2814 |
9.3137 |
19880601 |
20031231 |
24.06 |
|
Stuttgart-Echterdingen |
Baden-W<fc>rttemberg |
4931 |
371 |
48.6883 |
9.2235 |
19530101 |
20211231 |
25.61 |
Identify suitable weather stations for patching gaps
There are weather stations in the DWD list that are mostly complete I think, but for the sake of me trying to create a list (and we can never be sure that they indeed contain complete data with just the range of years shown), I will try to use 3 stations: 1, 3, and 6 from the weather station list.
Download weather data for promising stations, convert them to
chillR format and compile them in a list
##attempt 3
require(reshape2)
require(kableExtra)
require(ggplot2)
library(chillR)
station_list_dwd<-handle_dwd(action="list_stations",location=c(9.40,48.49),time_interval=c(19900101,20201231))
positions_in_station_list<-c(1,3,6)
patch_weather<-list()
for(i in 1:length(positions_in_station_list))
{patch_weather[[i]]<-handle_dwd(handle_dwd(action="download_weather",
location=station_list_dwd$Station_ID[positions_in_station_list[i]],time_interval=c(19900101,20211231)))[[1]]
names(patch_weather)[i]<-station_list_dwd$Station_name[positions_in_station_list[i]]
}
save_temperature_scenarios(patch_weather,"data/", "patch_weather")Use the patch_daily_temperatures function to fill
gaps
Instead of Chapter 30, this exercise is the last bit I worked on
because I had been wracking my head on how to carry this out when the
patch_daily_temps seemed to only work perfectly fine with
‘gsod’ datafiles. So many hiccups because I could not also clean the
‘gsod’ raw files with handle_gsod for some reason. It
always returns an empty file. That was the reason I had to use ‘dwd’
datafiles and move the location to Swabian Orchard paradise.
After so many attempts, I was able to figure out that the ‘old and
deprecated’ function patch_daily_temperatures would be able
to handle ‘dwd’ datafiles well– at least in my case. Please remember
that the patched data here is just for the sake of doing this exercise.
I will be using the complete authentic data from dwd for my location for
the next exercises. So, here’s the product of my many hours of sweat and
tears:
require(reshape2)
require(kableExtra)
require(ggplot2)
library(chillR)
##recall previous list()
station_list_dwd<-handle_dwd(action="list_stations",location=c(9.40,48.49),time_interval=c(19900101,20201231))
positions_in_station_list<-c(1,3,6)
patch_weather<-list()
for(i in 1:length(positions_in_station_list))
{patch_weather[[i]]<-handle_dwd(handle_dwd(action="download_weather",
location=station_list_dwd$Station_ID[positions_in_station_list[i]],time_interval=c(19900101,20211231)))[[1]]
names(patch_weather)[i]<-station_list_dwd$Station_name[positions_in_station_list[i]]
}
##assign the dataset with gaps to object
MApfelstetten<-read_tab("data/MApfelstetten_chillR_weather_gaps.csv")
#first patching, I used the old function 'patch_daily_temperatures'
patched<-patch_daily_temperatures(weather = MApfelstetten,
patch_weather = patch_weather)
##to summarize statistics:
# > patched$statistics
#$`Lenningen-Schopfloch`
# mean_bias stdev_bias filled gaps_remain
#Tmin -1.709 2.032 80 1027
#Tmax -0.085 1.283 80 1027
#$`Reutlingen-Betzingen`
# mean_bias stdev_bias filled gaps_remain
#Tmin -2.062 1.664 0 1027
#Tmax -2.999 1.501 0 1027
#$`Stuttgart-Echterdingen`
# mean_bias stdev_bias filled gaps_remain
#Tmin -1.950 1.595 1027 0
#Tmax -2.352 1.383 1027 0
##thus no need to limit mean bias to 1 and stdevbias to 2
write.csv(patched$weather, "data/MApfelstetten_patched.csv", row.names = FALSE)Investigate the results - have all gaps been filled?
library(chillR)
patched<-read.csv("data/MApfelstetten_patched.csv")
post_patch_stats<-fix_weather(patched)$QC
write.csv(post_patch_stats, "data/MApfelstetten_post_patchstats.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Missing_Tmin |
Missing_Tmax |
Incomplete_days |
Perc_complete |
|---|---|---|---|---|---|---|---|
|
1989/1990 |
1990 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1990/1991 |
1991 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1991/1992 |
1992 |
366 |
366 |
0 |
0 |
0 |
100 |
|
1992/1993 |
1993 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1993/1994 |
1994 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1994/1995 |
1995 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1995/1996 |
1996 |
366 |
366 |
0 |
0 |
0 |
100 |
|
1996/1997 |
1997 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1997/1998 |
1998 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1998/1999 |
1999 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1999/2000 |
2000 |
366 |
366 |
0 |
0 |
0 |
100 |
|
2000/2001 |
2001 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2001/2002 |
2002 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2002/2003 |
2003 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2003/2004 |
2004 |
366 |
366 |
0 |
0 |
0 |
100 |
|
2004/2005 |
2005 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2005/2006 |
2006 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2006/2007 |
2007 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2007/2008 |
2008 |
366 |
366 |
0 |
0 |
0 |
100 |
|
2008/2009 |
2009 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2009/2010 |
2010 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2010/2011 |
2011 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2011/2012 |
2012 |
366 |
366 |
0 |
0 |
0 |
100 |
|
2012/2013 |
2013 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2013/2014 |
2014 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2014/2015 |
2015 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2015/2016 |
2016 |
366 |
366 |
0 |
0 |
0 |
100 |
|
2016/2017 |
2017 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2017/2018 |
2018 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2018/2019 |
2019 |
365 |
365 |
0 |
0 |
0 |
100 |
|
2019/2020 |
2020 |
366 |
366 |
0 |
0 |
0 |
100 |
If necessary, repeat until you have a dataset you can work with in further analyses
Since everything has been filled, no need to repeat procedure. What I wish to repeat though is that I will not use the produced dataset for further exercises.
Exercise:
For the location you chose for your earlier analyses, use
chillR’s weather generator to produce 100 years of
synthetic temperature data.
I will be using Swabian Orchard Paradise (48.4931° N, 9.3998° E) for this portion.
library(chillR)
library(ggplot2)
MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
Temp<-temperature_generation(MApfelstetten_weather,
years=c(1998,2005),
sim_years = c(2001,2100))
Temperatures<-cbind(MApfelstetten_weather[
which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
Temperatures<-rbind(Temperatures,
cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
Data_source="simulated"))
Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
# chill_comparison_full_seasons<-read.csv("data/chill_comparison_full_seasons.csv")
ggplot(data=Temperatures, aes(Date,Tmin)) +
geom_smooth(aes(colour = factor(Year))) +
facet_wrap(vars(Data_source)) +
theme_bw(base_size = 20) +
theme(legend.position = "none") +
scale_x_date(date_labels = "%b")
ggplot(data=Temperatures, aes(Date,Tmax)) +
geom_smooth(aes(colour = factor(Year))) +
facet_wrap(vars(Data_source)) +
theme_bw(base_size = 20) +
theme(legend.position = "none") +
scale_x_date(date_labels = "%b")
Calculate winter chill (in Chill Portions) for your synthetic weather, and illustrate your results as histograms and cumulative distributions.
library(chillR)
library(ggplot2)
MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
Temp<-temperature_generation(MApfelstetten_weather,
years=c(1998,2005),
sim_years = c(2001,2100))
Temperatures<-cbind(MApfelstetten_weather[
which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
Temperatures<-rbind(Temperatures,
cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
Data_source="simulated"))
Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
chill_observed<-chilling(
stack_hourly_temps(
Temperatures[which(Temperatures$Data_source=="observed"),],
latitude = 48.49),
Start_JDay = 305,
End_JDay = 59)
chill_simulated<-chilling(
stack_hourly_temps(
Temperatures[which(Temperatures$Data_source=="simulated"),],
latitude = 48.49),
Start_JDay = 305,
End_JDay = 59)
chill_comparison<-cbind(chill_observed ,Data_source="observed")
chill_comparison<-rbind(chill_comparison,
cbind(chill_simulated ,Data_source="simulated"))
chill_comparison_full_seasons<-chill_comparison[
which(chill_comparison$Perc_complete==100),]
ggplot(chill_comparison_full_seasons, aes(x=Chill_portions)) +
geom_histogram(binwidth=1,aes(fill = factor(Data_source))) +
theme_bw(base_size = 20) +
labs(fill = "Data source") +
xlab("Chill accumulation (Chill Portions)") +
ylab("Frequency")
chill_simulations<-chill_comparison_full_seasons[
which(chill_comparison_full_seasons$Data_source=="simulated"),]
ggplot(chill_simulations, aes(x=Chill_portions)) +
stat_ecdf(geom = "step",lwd=1.5,col="blue") +
ylab("Cumulative probability") +
xlab("Chill accumulation (in Chill Portions)") +
theme_bw(base_size = 20)
write.csv(chill_comparison_full_seasons,"data/chill_comparison_full_seasons.csv", row.names = FALSE)
|
Season |
End_year |
Season_days |
Data_days |
Perc_complete |
Chilling_Hours |
Utah_Model |
Chill_portions |
GDH |
Data_source |
|---|---|---|---|---|---|---|---|---|---|
|
1998/1999 |
1999 |
120 |
120 |
100 |
1016 |
700.5 |
57.31928 |
497.6681 |
observed |
|
1999/2000 |
2000 |
120 |
120 |
100 |
1209 |
810.0 |
62.26994 |
632.8378 |
observed |
|
2000/2001 |
2001 |
121 |
121 |
100 |
1621 |
1262.5 |
74.80941 |
1027.6474 |
observed |
|
2001/2002 |
2002 |
120 |
120 |
100 |
1220 |
992.5 |
66.90289 |
855.7386 |
observed |
|
2002/2003 |
2003 |
120 |
120 |
100 |
1301 |
1055.0 |
66.67006 |
702.7373 |
observed |
|
2003/2004 |
2004 |
120 |
120 |
100 |
1211 |
845.0 |
65.66647 |
1081.6018 |
observed |
##amount of chill that is exceeded in 90% of all years
quantile(chill_simulations$Chill_portions, 0.1)
# 10%
#59.0096
##chill at 50% confidence interval (25th and 75th percentile)
quantile(chill_simulations$Chill_portions, c(0.25,0.75))
# 25% 75%
#62.81098 71.61968 Produce similar plots for the number of freezing hours (<0°C) in April (or October, if your site is in the Southern Hemisphere) for your location of interest.
#attempt 3, freezing() function created based on chillR's chilling function
library(chillR)
library(ggplot2)
##synthetic weather generation
MApfelstetten_weather<-read.csv("data/MApfelstetten_weather.csv")
Temp<-temperature_generation(MApfelstetten_weather,
years=c(1998,2005),
sim_years = c(2001,2100))
Temperatures<-cbind(MApfelstetten_weather[
which(MApfelstetten_weather$Year %in% 1998:2005),] ,Data_source="observed")
Temperatures<-rbind(Temperatures,
cbind(Temp[[1]][,c("Year","Month","Day","Tmin","Tmax")],
Data_source="simulated"))
Temperatures[,"Date"]<-as.Date(ISOdate(2000, Temperatures$Month, Temperatures$Day))
##creating function calculating freezing hours based on chillR function chilling()
freezing<-function (hourtemps = NULL, Start_JDay = 1, End_JDay = 366, THourly = NULL,
misstolerance = 50)
{if (is.null(hourtemps) & !is.null(THourly))
hourtemps <- THourly
if ((length(names(hourtemps)) == 2) & ("hourtemps" %in% names(hourtemps))) {
QC <- hourtemps$QC
hourtemps <- hourtemps$hourtemps
}
else QC <- NULL
if (Start_JDay < End_JDay) {
hourtemps[which(hourtemps$JDay >= Start_JDay & hourtemps$JDay <=
End_JDay), "sea"] <- hourtemps[which(hourtemps$JDay >=
Start_JDay & hourtemps$JDay <= End_JDay), "Year"]
}
else {
hourtemps[which(hourtemps$JDay >= Start_JDay), "sea"] <- hourtemps[which(hourtemps$JDay >=
Start_JDay), "Year"] + 1
hourtemps[which(hourtemps$JDay <= End_JDay), "sea"] <- hourtemps[which(hourtemps$JDay <=
End_JDay), "Year"]
}
if (Start_JDay < End_JDay) {
relevant_days <- Start_JDay:End_JDay
}
else {
relevant_days <- c(Start_JDay:366, 1:End_JDay)
}
normal_lines <- which(!(hourtemps$JDay == Start_JDay & hourtemps$Hour ==
0))
normal_lines <- normal_lines[which(normal_lines > 1)]
hourtemps <- hourtemps[which(!is.na(hourtemps[, "Temp"])),
]
##range for freezing hour counting, creating columns
FH_range <- which(hourtemps$Temp <
0)
hourtemps[, "FH_weights"] <- 0
hourtemps[FH_range, "FH_weights"] <- 1
hourtemps[, "FH"] <- 0
seasons <- as.numeric(unique(hourtemps$sea))
seasons <- seasons[!is.na(seasons)]
##storing seasons
chillout <- data.frame(Season = paste(seasons - 1, "/", seasons,
sep = ""), End_year = seasons)
if (End_JDay >= Start_JDay)
dt <- sapply(seasons, function(x) difftime(ISOdate(x -
1, 12, 31) + End_JDay * 86400, ISOdate(x - 1, 12,
31) + Start_JDay * 86400))
if (End_JDay < Start_JDay)
dt <- sapply(seasons, function(x) difftime(ISOdate(x -
1, 12, 31) + End_JDay * 86400, ISOdate(x - 2, 12,
31) + Start_JDay * 86400))
chillout[, "Season_days"] <- dt + 1
chillout[, "Data_days"] <- sapply(seasons, function(x) length(which(hourtemps$sea ==
x))/24)
for (sea in seasons) {
if ("no_Tmin" %in% names(hourtemps) & "no_Tmax" %in%
names(hourtemps))
chillout[which(chillout$End_year == sea), "Interpolated_days"] <- length(which(hourtemps[which(hourtemps$sea ==
sea), "no_Tmin"] | hourtemps[which(hourtemps$sea ==
sea), "no_Tmax"]))/24
chillout[, "Perc_complete"] <- NA
##creation of column for summed freezing hours
chillout[which(chillout$End_year == sea), "Freezing_Hours"] <- sum(hourtemps[which(hourtemps$sea ==
sea), "FH_weights"])
}
if ("no_Tmin" %in% names(hourtemps) & "no_Tmax" %in% names(hourtemps))
chillout[, "Perc_complete"] <- (chillout[, "Data_days"] -
chillout[, "Interpolated_days"])/chillout[, "Season_days"] *
100
else chillout[, "Perc_complete"] <- chillout[, "Data_days"]/chillout[,
"Season_days"] * 100
chillout <- chillout[which(chillout$Perc_complete >= 100 -
misstolerance), ]
return(chillout)
}
freezing_observed<-freezing(
stack_hourly_temps(
Temperatures[which(Temperatures$Data_source=="observed"),],
latitude = 48.49),
Start_JDay = 91,
End_JDay = 120)
freezing_simulated<-freezing(
stack_hourly_temps(
Temperatures[which(Temperatures$Data_source=="simulated"),],
latitude = 48.49),
Start_JDay = 91,
End_JDay = 120)
freeze_comparison<-cbind(freezing_observed ,Data_source="observed")
freeze_comparison<-rbind(freeze_comparison,
cbind(freezing_simulated ,Data_source="simulated"))
freeze_comparison_full_seasons<-freeze_comparison[
which(freeze_comparison$Perc_complete==100),]
write.csv(freeze_comparison_full_seasons, "data/freeze_comparison_full_seasons.csv", row.names = FALSE)
##generating the Cumulative Freezing (sum of Freezing hours) vs Frequency plot
ggplot(freeze_comparison_full_seasons, aes(x=Freezing_Hours)) +
geom_histogram(binwidth=1,aes(fill = factor(Data_source))) +
theme_bw(base_size = 20) +
labs(fill = "Data source") +
xlab("Cumulative Freezing Hours") +
ylab("Frequency")
##generating the Cumulative Freezing (sum of Freezing hours) vs Cumulative Probability plot
freeze_simulations<-freeze_comparison_full_seasons[
which(freeze_comparison_full_seasons$Data_source=="simulated"),]
ggplot(freeze_simulations, aes(x=Freezing_Hours)) +
stat_ecdf(geom = "step",lwd=1.5,col="blue") +
ylab("Cumulative probability") +
xlab("Cumulative Freezing Hours") +
theme_bw(base_size = 20)
|
Season |
End_year |
Season_days |
Data_days |
Perc_complete |
Freezing_Hours |
Data_source |
|---|---|---|---|---|---|---|
|
1997/1998 |
1998 |
30 |
30 |
100 |
42 |
observed |
|
1998/1999 |
1999 |
30 |
30 |
100 |
48 |
observed |
|
1999/2000 |
2000 |
30 |
30 |
100 |
42 |
observed |
|
2000/2001 |
2001 |
30 |
30 |
100 |
117 |
observed |
|
2001/2002 |
2002 |
30 |
30 |
100 |
48 |
observed |
|
2002/2003 |
2003 |
30 |
30 |
100 |
142 |
observed |
##amount of chill that is exceeded in 90% of all years
quantile(freeze_simulations$Freezing_Hours, 0.1)
#10%
#7.9
##chill at 50% confidence interval (25th and 75th percentile)
quantile(freeze_simulations$Freezing_Hours, c(0.25,0.75))
#25% 75%
#23.00 99.25 Exercise:
For the location you chose for previous exercises, produce historic temperature scenarios representing several years of the historic record (your choice).
After looking at the mean annual temperatures in Germany in a paper by Glaser et. al (2017), I saw that the temperatures only started to rise steeply during 1980s (same as what we had in the samples), so I decided to take a look at those same years too. Table from the paper shown below:
require(chillR)
require(ggplot2)
require(reshape2)
require(kableExtra)
#download station list for time interval of interest
station_list<-handle_dwd(action="list_stations",
location=c(9.40,48.49),
time_interval = c(19730101,20191231))
station_list
# download weather data from M<fc>nsingen-Apfelstetten (2nd in list) since it is still in the list and I've dealt with it before
MApfelstetten_weather_raw<-handle_dwd(action="download_weather",
location=station_list$Station_ID[2],
time_interval = c(19730101,20191231),
station_list = station_list)
# convert weather data to chillR format
MApfelstetten_weather<-handle_dwd(MApfelstetten_weather_raw)
# check record for missing data
MApfelstetten_historic<-fix_weather(MApfelstetten_weather[[1]])$QC
write.csv(MApfelstetten_historic, "data/MApfelstetten_historic.csv", row.names = FALSE)
##there is no gap for the data set, proceeding to generating scenarios for 1980s
MApfelstetten_temps<-MApfelstetten_weather[[1]]
scenario_1980<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1980)
scenario_1980$'1980'$data
scenario_1980
temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = scenario_1980)
scenario_1996<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1996)
relative_scenario<-temperature_scenario_baseline_adjustment(
baseline=scenario_1996,
temperature_scenario = scenario_1980)
temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = relative_scenario)
all_past_scenarios<-temperature_scenario_from_records(
weather=MApfelstetten_temps,
year=c(1980,1990,2000,2010))
adjusted_scenarios<-temperature_scenario_baseline_adjustment(
baseline=scenario_1996,
temperature_scenario = all_past_scenarios)
all_past_scenario_temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = adjusted_scenarios)
all_past_scenario_temps
|
Season |
End_year |
Season_days |
Data_days |
Missing_Tmin |
Missing_Tmax |
Incomplete_days |
Perc_complete |
|---|---|---|---|---|---|---|---|
|
1972/1973 |
1973 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1973/1974 |
1974 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1974/1975 |
1975 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1975/1976 |
1976 |
366 |
366 |
0 |
0 |
0 |
100 |
|
1976/1977 |
1977 |
365 |
365 |
0 |
0 |
0 |
100 |
|
1977/1978 |
1978 |
365 |
365 |
0 |
0 |
0 |
100 |
Produce chill distributions for these scenarios and plot them.
require(chillR)
require(ggplot2)
require(reshape2)
require(kableExtra)
##copying and pasting the code from previous step
#download station list for time interval of interest
station_list<-handle_dwd(action="list_stations",
location=c(9.40,48.49),
time_interval = c(19730101,20191231))
station_list
# download weather data from M<fc>nsingen-Apfelstetten (2nd in list) since it is still in the list and I've dealt with it before
MApfelstetten_weather_raw<-handle_dwd(action="download_weather",
location=station_list$Station_ID[2],
time_interval = c(19730101,20191231),
station_list = station_list)
# convert weather data to chillR format
MApfelstetten_weather<-handle_dwd(MApfelstetten_weather_raw)
# check record for missing data
MApfelstetten_historic<-fix_weather(MApfelstetten_weather[[1]])$QC
write.csv(MApfelstetten_historic, "data/MApfelstetten_historic.csv", row.names = FALSE)
##there is no gap for the data set, proceeding to generating scenarios for 1980s
MApfelstetten_temps<-MApfelstetten_weather[[1]]
write.csv(MApfelstetten_temps, "data/MApfelstetten_temps.csv", row.names = FALSE)
scenario_1980<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1980)
scenario_1980$'1980'$data
scenario_1980
temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = scenario_1980)
scenario_1996<-temperature_scenario_from_records(weather=MApfelstetten_temps,year=1996)
relative_scenario<-temperature_scenario_baseline_adjustment(
baseline=scenario_1996,
temperature_scenario = scenario_1980)
temps_1980<-temperature_generation(weather=MApfelstetten_temps, years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = relative_scenario)
all_past_scenarios<-temperature_scenario_from_records(
weather=MApfelstetten_temps,
year=c(1980,1990,2000,2010))
adjusted_scenarios<-temperature_scenario_baseline_adjustment(
baseline=scenario_1996,
temperature_scenario = all_past_scenarios)
all_past_scenario_temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2100),
temperature_scenario = adjusted_scenarios)
all_past_scenario_temps
#chill accumulation calculation
chill_hist_scenario_list<-tempResponse_daily_list(all_past_scenario_temps,latitude=48.49,Start_JDay = 305, End_JDay = 59)
chill_hist_scenario_list
#converting list to data.frame
scenarios<-names(chill_hist_scenario_list)[1:4]
all_scenarios<-chill_hist_scenario_list[[scenarios[1]]]
all_scenarios[,"scenario"]<-as.numeric(scenarios[1])
for (sc in scenarios[2:4])
all_scenarios<-rbind(all_scenarios,
cbind(chill_hist_scenario_list[[sc]],scenario=as.numeric(sc)))
all_scenarios<-all_scenarios[which(all_scenarios$Perc_complete==100),]
#compute the actual 'observed' chill for comparison
actual_chill<-tempResponse_daily_list(MApfelstetten_temps,latitude=48.49,
Start_JDay = 305,
End_JDay = 59)[[1]]
actual_chill<-actual_chill[which(actual_chill$Perc_complete==100),]
actual_chill
ggplot(data=all_scenarios,aes(scenario,Chill_Portions,
fill=factor(scenario))) +
geom_violin() +
ylab("Chill accumulation (Chill Portions)") +
xlab("Scenario year") +
theme_bw(base_size=12) +
ylim(c(0,90)) +
geom_point(data=actual_chill,
aes(End_year,Chill_Portions,fill="black"),
col="blue",show.legend = FALSE) +
scale_fill_discrete(name="Scenario",
breaks = unique(all_scenarios$scenario))
ggplot(data=all_scenarios,aes(scenario,GDH,
fill=factor(scenario))) +
geom_violin() +
ylab("Heat accumulation (Growing Degree Hours)") +
xlab("Scenario year") +
theme_bw(base_size=12) +
ylim(c(0,7000)) +
geom_point(data=actual_chill,
aes(End_year,GDH,fill="black"),
col="blue",show.legend = FALSE) +
scale_fill_discrete(name="Scenario",
breaks = unique(all_scenarios$scenario))
?runn_mean
temperature_means<-data.frame(Year=min(MApfelstetten_temps$Year):max(MApfelstetten_temps$Year),
Tmin=aggregate(MApfelstetten_temps$Tmin,FUN="mean",
by=list(MApfelstetten_temps$Year))[,2],
Tmax=aggregate(MApfelstetten_temps$Tmax,FUN="mean",
by=list(MApfelstetten_temps$Year))[,2])
temperature_means[,"runn_mean_Tmin"]<-runn_mean(temperature_means$Tmin,15)
temperature_means[,"runn_mean_Tmax"]<-runn_mean(temperature_means$Tmax,15)
Tmin_regression<-lm(Tmin~Year, temperature_means)
temperature_means[,"regression_Tmin"]<-Tmin_regression$coefficients[1]+
Tmin_regression$coefficients[2]*temperature_means$Year
Tmax_regression<-lm(Tmax~Year, temperature_means)
temperature_means[,"regression_Tmax"]<-Tmax_regression$coefficients[1]+
Tmax_regression$coefficients[2]*temperature_means$Year
ggplot(temperature_means,aes(Year, Tmin)) +
geom_point() +
geom_line(data=temperature_means,aes(Year, runn_mean_Tmin),lwd=2,col="blue") +
geom_line(data=temperature_means,aes(Year, regression_Tmin),lwd=2,col="red") +
theme_bw(base_size=12) +
ylab("Mean monthly minimum temperature (°C)")
ggplot(temperature_means,aes(Year, Tmax)) +
geom_point() +
geom_line(data=temperature_means,aes(Year, runn_mean_Tmax),lwd=2,col="blue") +
geom_line(data=temperature_means,aes(Year, regression_Tmax),lwd=2,col="red") +
theme_bw(base_size=12) +
ylab("Mean monthly maximum temperature (°C)")
Exercise:
Analyze the historic and future impact of climate change on three agroclimatic metrics of your choice, for the location you’ve chosen for your earlier analyses.
install.packages("ggplot2")
require(chillR)
require(ggplot2)
require(reshape2)
require(kableExtra)
MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
RCPs<-c("rcp45,rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{start_year <- Time-15
end_year <- Time+15
clim_scen <-getClimateWizardData(
c(longitude = 9.399,latitude = 48.493),
RCP,
start_year,
end_year,
temperature_generation_scenarios = TRUE,
baseline =c(1975, 2005),
metric = "monthly_min_max_temps",
GCMs = "all")
save_temperature_scenarios(clim_scen,
"data/ClimateWizard",
paste0("MApfelstetten_futures_",Time,"_",RCP))}
scenario_1990<-temperature_scenario_from_records(MApfelstetten_temps,1990)
scenario_1996<-temperature_scenario_from_records(MApfelstetten_temps,1996)
adjustment_scenario<-temperature_scenario_baseline_adjustment(scenario_1996,scenario_1990)
RCPs<-c("rcp85")
Times<-c(2090)
for(RCP in RCPs)
for(Time in Times)
{
clim_scen<-load_ClimateWizard_scenarios(
"data/climateWizard",
paste0("MApfelstetten_futures_",2090,"_","rcp85"))
clim_scen_adjusted<-
temperature_scenario_baseline_adjustment(baseline_temperature_scenario = adjustment_scenario, temperature_scenario = clim_scen)
Temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2101),
temperature_scenario = clim_scen_adjusted)
save_temperature_scenarios(
Temps,
"data/Weather",
paste0("MApfelstetten_",2090,"_","rcp85"))
}
scenario_1990<-temperature_scenario_from_records(MApfelstetten_temps,1990)
scenario_1996<-temperature_scenario_from_records(MApfelstetten_temps,1996)
adjustment_scenario<-temperature_scenario_baseline_adjustment(scenario_1996,scenario_1990)
clim_scen<-load_ClimateWizard_scenarios(
"data/climateWizard",
paste0("MApfelstetten_futures_",2090,"_","rcp45"))
clim_scen
clim_scen_adjusted<-
temperature_scenario_baseline_adjustment(baseline_temperature_scenario = adjustment_scenario, temperature_scenario = clim_scen)
Temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2101),
temperature_scenario = clim_scen_adjusted)
save_temperature_scenarios(
Temps,
"data/Weather",
paste0("MApfelstetten_",2090,"_","rcp45"))
RCPs<-c("rcp85")
Times<-c(2085)
for(RCP in RCPs)
for(Time in Times)
{
clim_scen<-load_ClimateWizard_scenarios(
"data/climateWizard",
paste0("MApfelstetten_futures_",Time,"_",RCP))
clim_scen_adjusted<-
temperature_scenario_baseline_adjustment(
baseline_temperature_scenario=adjustment_scenario,
temperature_scenario=clim_scen)
Temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2101),
temperature_scenario = clim_scen_adjusted)
save_temperature_scenarios(
Temps,
"data/Weather",
paste0("MApfelstetten_",Time,"_",RCP))
}
##
all_past_scenarios<-temperature_scenario_from_records(
weather=MApfelstetten_temps,
year=c(1980,1990,2000,2010))
adjusted_scenarios<-temperature_scenario_baseline_adjustment(
baseline=scenario_1996,
temperature_scenario = all_past_scenarios)
all_past_scenario_temps<-temperature_generation(
weather=MApfelstetten_temps,
years=c(1973,2019),
sim_years=c(2001,2101),
temperature_scenario = adjusted_scenarios)
save_temperature_scenarios(
all_past_scenario_temps,
"data/Weather",
"MApfelstetten_historic")
frost_model<-function(x) step_model(x,data.frame(
lower=c(-1000,0),
upper=c(0,1000),
weight=c(1,0)))
models<-list(Chill_CP=Dynamic_Model,
Heat_GDH=GDH,
Frost_H=frost_model)
Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
chill_past_scenarios<-tempResponse_daily_list(
Temps,
latitude=48.49,
Start_JDay = 305,
End_JDay = 59,
models=models,
misstolerance = 10)
chill_observed<-tempResponse_daily_list(
MApfelstetten_temps,
latitude=48.49,
Start_JDay = 305,
End_JDay = 59,
models=models,
misstolerance = 10)
save_temperature_scenarios(chill_past_scenarios,
"data/chill",
"MApfelstetten_historic")
save_temperature_scenarios(chill_observed,
"data/chill",
"MApfelstetten_observed")
chill_past_scenarios<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_historic")
chill_observed<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_observed")
chills <-make_climate_scenario(
chill_past_scenarios,
caption = "Historic",
historic_data = chill_observed,
time_series = TRUE)
plot_climate_scenarios(
climate_scenario_list=chills,
metric="Chill_CP",
metric_label="Chill (Chill Portions)")
RCPs<-c("rcp45", "rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{
Temps<-load_temperature_scenarios(
"data/Weather",
paste0("MApfelstetten_",Time,"_",RCP))
chill<-tempResponse_daily_list(
Temps,
latitude=48.49,
Start_JDay = 305,
End_JDay = 59,
models=models,
misstolerance = 10)
save_temperature_scenarios(
chill,
"data/chill",
paste0("MApfelstetten_",Time,"_",RCP))
}
for(RCP in RCPs)
for(Time in Times)
{
chill<-load_temperature_scenarios(
"data/chill",
paste0("MApfelstetten_",Time,"_",RCP))
if(RCP=="rcp45") RCPcaption <- "RCP4.5"
if(RCP=="rcp85") RCPcaption <- "RCP8.5"
if(Time=="2050") Time_caption <- "2050"
if(Time=="2090") Time_caption <- "2090"
chills <-make_climate_scenario(
chill,
caption =c(RCPcaption, Time_caption),
add_to = chills)
}
#plotting
plot_climate_scenarios(
climate_scenario_list=chills,
metric="Chill_CP",
metric_label="Chill (Chill Portions)",
texcex=1.5)
plot_climate_scenarios(
climate_scenario_list=chills,
metric="Frost_H",
metric_label="Frost hours",
texcex=1.5)
Chill Portions-
the Swabian Orchard Paradise will not have to worry about not being able
to grow their trees up to 2090, it seems, even with RCP8.5
projection.
Heat (GDH)- per the projection, it will become warmer, thus can give opportunity for farmers to venture into more diverse tree species that require higher GDH in the future.
Frost hours- hours of frost is projected to decline, which is also good because less frost damage to flowers during critical times.
In general, these plots come in agreement with what we are expecting: that the location (which is in Germany) will be more fit for production- this was in line with what lecturers from other courses were saying. Global North will be benefiting more from climate change than the Global South.
Exercise: Produce similar plots for the weather station you selected for earlier exercises.
##install.packages("ggpmisc")
library(chillR)
library(ggplot2)
library(reshape2)
library(kableExtra)
library(ggpmisc)
library(patchwork)
chill_past_scenarios<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_historic")
chill_observed<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_observed")
chills <-make_climate_scenario(
chill_past_scenarios,
caption = "Historic",
historic_data = chill_observed,
time_series = TRUE)
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{ chill<-load_temperature_scenarios(
"data/chill",
paste0("MApfelstetten_",Time,"_",RCP))
if(RCP=="rcp45") RCPcaption <- "RCP4.5"
if(RCP=="rcp85") RCPcaption <- "RCP8.5"
if(Time=="2050") Time_caption <- "2050"
if(Time=="2090") Time_caption <- "2090"
chills <-make_climate_scenario(
chill,
caption =c(RCPcaption, Time_caption),
add_to = chills)
}
plot_climate_scenarios(
climate_scenario_list=chills,
metric="Frost_H",
metric_label="Frost hours",
texcex=1.5)
# We'll first process the past scenarios (element 1 of the chills list).
# Within the data element, we have a list of multiple data.frames for
# the various past scenarios.
# Using a 'for' loop, we cycle through all these data.frames.
for(nam in names(chills[[1]]$data))
{
# Extract the data frame.
ch<-chills[[1]]$data[[nam]]
# Add columns for the new information we have to add and fill them.
ch[,"GCM"]<-"none"
ch[,"RCP"]<-"none"
ch[,"Year"]<-as.numeric(nam)
# Now check if this is the first time we've gone through this loop.
# If this is the first time, the ch data.frame becomes the output
# object (past_simulated).
# If it is not the first time ('else'), we add the current data.frame
# to the 'past_simulated' object
if(nam==names(chills[[1]]$data)[1])
past_simulated<-ch else
past_simulated<-rbind(past_simulated,ch)
}
# We add another column called 'Scenario' and label all rows as 'Historic'
past_simulated["Scenario"] <- "Historic"
kable(past_simulated[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
# We'll want to add the historic observation too, so let's simplify the
# pointer to this information for easier use later
past_observed <- chills[[1]][["historic_data"]]
kable(past_observed[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
# Extract future data
for(i in 2:length(chills))
for(nam in names(chills[[i]]$data))
{ch<-chills[[i]]$data[[nam]]
ch[,"GCM"]<-nam
ch[,"RCP"]<-chills[[i]]$caption[1]
ch[,"Year"]<-chills[[i]]$caption[2]
if(i==2&nam==names(chills[[i]]$data)[1])
future_data<-ch else
future_data<-rbind(future_data,ch)
}
kable(future_data[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
# Extract the model names
Models <- unique(future_data$GCM)
metric<-"Heat_GDH"
axis_label<-"Heat (in GDH)"
# get extreme values for the axis scale
rng<-range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
rng
past_plot<-ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",metric,group="Year"),
fill="skyblue")
past_plot
past_plot<-past_plot +
scale_y_continuous(limits = c(0, round(rng[2] + rng[2]/10))) +
labs(x = "Year", y = axis_label)
past_plot
past_plot<-past_plot +
facet_grid(~ Scenario) +
theme_bw(base_size = 15)
past_plot
past_plot<-past_plot +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1))
past_plot
# add historic data
past_plot <- past_plot +
geom_point(data = past_observed,
aes_string("End_year",metric),
col="blue")
past_plot
y<-2050
future_2050<-ggplot(data= future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM", metric, fill="GCM"))
future_2050
future_2050<-future_2050 +
facet_wrap(vars(RCP)) +
scale_x_discrete(labels = NULL, expand = expansion(add = 1))
##stopped here
future_2050<-future_2050 +
scale_y_continuous(limits = c(0, round(1.1*rng[2]))) +
geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5)
future_2050
future_2050<-future_2050 +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
face = "bold"))
future_2050
future_plot_list<-list()
for(y in c(2050,2090))
{
future_plot_list[[which(y == c(2050,2090))]] <-
ggplot(data= future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
facet_wrap(vars(RCP)) +
scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center", npcy = "top", label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(
hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
future_plot_list
library(patchwork)
both_plots<-past_plot + future_plot_list
both_plots
plot <- both_plots +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
plot<-plot & theme(legend.position = "bottom",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
metric<-"Chill_CP"
axis_label<-"Chill (in CP)"
# get extreme values for the axis scale
rng<-range(past_observed[[metric]],past_simulated[[metric]],future_data[[metric]])
past_plot<-ggplot() +
geom_boxplot(data = past_simulated, aes_string("as.numeric(Year)",metric,group="Year"),fill="skyblue") +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1)) +
geom_point(data = past_observed,aes_string("End_year",metric),col="blue")
future_plot_list<-list()
for(y in c(2050,2090))
{
future_plot_list[[which(y == c(2050,2090))]] <- ggplot(data= future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
facet_wrap(vars(RCP)) +
scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
face = "bold"))
}
plot<- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
metric<-"Frost_H"
axis_label<-"Frost duration (in hours)"
# get extreme values for the axis scale
rng<-range(past_observed[[metric]],past_simulated[[metric]],future_data[[metric]])
past_plot<-ggplot() +
geom_boxplot(data = past_simulated, aes_string("as.numeric(Year)",metric,group="Year"),fill="skyblue") +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1)) +
geom_point(data = past_observed,aes_string("End_year",metric),col="blue")
future_plot_list<-list()
for(y in c(2050,2090))
{
future_plot_list[[which(y == c(2050,2090))]] <- ggplot(data= future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
facet_wrap(vars(RCP)) +
scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center", npcy = "top", label = Year), size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5, vjust = -1, size = 15 * 1.05,
face = "bold"))
}
plot<- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
plot_scenarios_gg<-function(past_observed,
past_simulated,
future_data,
metric,
axis_label)
{
rng<-range(past_observed[[metric]],
past_simulated[[metric]],
future_data[[metric]])
past_plot<-ggplot() +
geom_boxplot(data = past_simulated,
aes_string("as.numeric(Year)",metric,group="Year"),
fill="skyblue") +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
labs(x = "Year", y = axis_label) +
facet_grid(~ Scenario) +
theme_bw(base_size = 15) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(angle=45, hjust=1)) +
geom_point(data = past_observed,
aes_string("End_year",metric),
col="blue")
future_plot_list<-list()
for(y in c(2050,2090))
{
future_plot_list[[which(y == c(2050,2090))]] <-
ggplot(data= future_data[which(future_data$Year==y),]) +
geom_boxplot(aes_string("GCM", metric, fill="GCM")) +
facet_wrap(vars(RCP)) +
scale_x_discrete(labels = NULL, expand = expansion(add = 1)) +
scale_y_continuous(limits = c(0, round(round(1.1*rng[2])))) +
geom_text_npc(aes(npcx = "center", npcy = "top", label = Year),
size = 5) +
theme_bw(base_size = 15) +
theme(axis.ticks.y = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
legend.position = "bottom",
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.background = element_rect(),
strip.background = element_blank(),
strip.text = element_text(face = "bold"),
legend.box.spacing = unit(0, "cm"),
plot.subtitle = element_text(hjust = 0.5,
vjust = -1,
size = 15 * 1.05,
face = "bold"))
}
plot<- (past_plot +
future_plot_list +
plot_layout(guides = "collect",
widths = c(1,rep(1.8,length(future_plot_list))))
) & theme(legend.position = "bottom",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
}
##plotting and saving plots (I manually created the "plots" folder)
plot_scenarios_gg(past_observed=past_observed,
past_simulated=past_simulated,
future_data=future_data,
metric="Heat_GDH",
axis_label="Heat (in Growing Degree Hours)")
ggsave("data/plots/Heat_GDH.png", width = 20,
height = 10, units = "cm",
dpi = 600)
plot_scenarios_gg(past_observed=past_observed,
past_simulated=past_simulated,
future_data=future_data,
metric="Chill_CP",
axis_label="Chill (in Chill Portions)")
ggsave("data/plots/Chill_CP.png", width = 20,
height = 10, units = "cm",
dpi = 600)
plot_scenarios_gg(past_observed=past_observed,
past_simulated=past_simulated,
future_data=future_data,
metric="Frost_H",
axis_label="Frost duration (in hours)")
ggsave("data/plots/Frost_H.png", width = 20,
height = 10, units = "cm",
dpi = 600)
Exercise:
Perform a similar analysis for the location you’ve chosen for your exercises.
I just want to share that I got super stuck in installing dormancyR because typing in the devtools library constantly gave an error that looked like this:
##Error: package or namespace load failed for ‘devtools’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
## namespace ‘htmltools’ 0.5.3 is already loaded, but >= 0.5.4 is required
I had to install and use the “remotes” package for it to work. Anyway, here’s the code:
library(chillR)
##install.packages("remotes")
##install.packages("colorRamps")
install.packages("gganimate")
library(remotes)
remotes::install_github("EduardoFernandezC/dormancyR")
library(dormancyR)
library(ggplot2)
library(reshape2)
library(colorRamps)
library(kableExtra)
library(patchwork)
library(gganimate)
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
hourly_models <- list(Chilling_units = chilling_units,
Low_chill = low_chill_model,
Modified_Utah = modified_utah_model,
North_Carolina = north_carolina_model,
Positive_Utah = positive_utah_model,
Chilling_Hours = Chilling_Hours,
Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model)
daily_models<-list(Rate_of_Chill = rate_of_chill,
Chill_Days = chill_days,
Exponential_Chill = exponential_chill,
Triangula_Chill_Haninnen = triangular_chill_1,
Triangular_Chill_Legave = triangular_chill_2)
metrics<-c(names(daily_models),names(hourly_models))
model_labels=c("Rate of Chill",
"Chill Days",
"Exponential Chill",
"Triangular Chill (Häninnen)",
"Triangular Chill (Legave)",
"Chilling Units",
"Low-Chill Chill Units",
"Modified Utah Chill Units",
"North Carolina Chill Units",
"Positive Utah Chill Units",
"Chilling Hours",
"Utah Chill Units",
"Chill Portions")
kable(data.frame(Metric=model_labels,'Function name'=metrics)) %>%
kable_styling("striped", position = "left",font_size = 10)
MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
Start_JDay<-305
End_JDay<-59
daily_models_past_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
hourly_models_past_scenarios<-tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
#past_scenarios<-daily_models_past_scenarios
past_scenarios<-lapply(
names(daily_models_past_scenarios),
function(x)
cbind(daily_models_past_scenarios[[x]],
hourly_models_past_scenarios[[x]][,names(hourly_models)]))
names(past_scenarios)<-names(daily_models_past_scenarios)
daily_models_observed<-tempResponse_daily(
MApfelstetten_temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_observed<-
daily_models_observed[which(daily_models_observed$Perc_complete>90),]
hourly_models_observed<-tempResponse_daily_list(
MApfelstetten_temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
past_observed<-cbind(
daily_models_observed,
hourly_models_observed[[1]][,names(hourly_models)])
save_temperature_scenarios(past_scenarios,
"data/chill",
"MApfelstetten_multichill_historic")
write.csv(past_observed,
"data/chill/MApfelstetten_multichill_observed.csv",
row.names=FALSE)
##same procedure for future data
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{
Temps<-load_temperature_scenarios(
"data/Weather",
paste0("MApfelstetten_",Time,"_",RCP))
daily_models_future_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_future_scenarios<-lapply(
daily_models_future_scenarios,
function(x) x[which(x$Perc_complete>90),])
hourly_models_future_scenarios<-
tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
future_scenarios<-daily_models_future_scenarios
future_scenarios<-lapply(
names(future_scenarios),
function(x)
cbind(future_scenarios[[x]],
hourly_models_future_scenarios[[x]][,names(hourly_models)]))
names(future_scenarios)<-names(daily_models_future_scenarios)
chill<-future_scenarios
save_temperature_scenarios(
chill,
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
}
##loading and producing chill scenarios
chill_past_scenarios<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_multichill_historic")
chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
chills <-make_climate_scenario(chill_past_scenarios,
caption = "Historic",
historic_data = chill_observed,
time_series = TRUE)
for(RCP in RCPs)
for(Time in Times)
{
chill<-load_temperature_scenarios(
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
if(RCP=="rcp45") RCPcaption <- "RCP4.5"
if(RCP=="rcp85") RCPcaption <- "RCP8.5"
if(Time=="2050") Time_caption <- "2050"
if(Time=="2090") Time_caption <- "2090"
chills <-make_climate_scenario(chill,
caption =c(RCPcaption,Time_caption),
add_to = chills)
}
##heat map for Safe Winter chill
for(i in 1:length(chills))
{ch<-chills[[i]]
if(ch$caption[1]=="Historic")
{GCMs<-rep("none",length(names(ch$data)))
RCPs<-rep("none",length(names(ch$data)))
Years<-as.numeric(ch$labels)
Scenario<-rep("Historic",length(names(ch$data)))} else
{GCMs<-names(ch$data)
RCPs<-rep(ch$caption[1],length(names(ch$data)))
Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
Scenario<-rep("Future",length(names(ch$data)))}
for(nam in names(ch$data))
{for(met in metrics)
{temp_res<-data.frame(Metric=met,
GCM=GCMs[which(nam==names(ch$data))],
RCP=RCPs[which(nam==names(ch$data))],
Year=Years[which(nam==names(ch$data))],
Result=quantile(ch$data[[nam]][,met],0.1),
Scenario=Scenario[which(nam==names(ch$data))])
if(i==1&nam==names(ch$data)[1]&met==metrics[1])
results<-temp_res else
results<-rbind(results,temp_res)
}}}
for(met in metrics)
results[which(results$Metric==met),"SWC"]<-
results[which(results$Metric==met),"Result"]/
results[which(results$Metric==met&results$Year==1980),"Result"]-1
rng = range(results$SWC)
p_future<-ggplot(results[which(!results$GCM=="none"),],
aes(GCM, y=factor(Metric, levels=metrics),
fill = SWC)) +
geom_tile()
p_future
p_future <-
p_future +
facet_grid(RCP ~ Year)
p_future
#make text smaller
p_future <-
p_future +
theme_bw(base_size = 15) +
theme(axis.text = element_text(size=8))
p_future
library(colorRamps)
p_future <-
p_future +
scale_fill_gradientn(colours=matlab.like(15),
labels = scales::percent,
limits=rng)
p_future
p_future <-
p_future +
theme(axis.text.x = element_text(angle = 75, hjust = 1, vjust = 1)) +
labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
scale_y_discrete(labels=model_labels) +
ylab("Chill metric")
p_future
#plot of past scenarios
p_past<-
ggplot(results[which(results$GCM=="none"),],
aes(Year, y=factor(Metric, levels=metrics),
fill = SWC)) +
geom_tile()
p_past<-
p_past +
theme_bw(base_size = 12) +
theme(axis.text = element_text(size=8))
p_past<-
p_past +
scale_fill_gradientn(colours=matlab.like(15),
labels = scales::percent,
limits=rng)
p_past<-
p_past +
scale_x_continuous(position = "top")
p_past<-
p_past +
labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
scale_y_discrete(labels=model_labels) +
ylab("Chill metric")
chill_comp_plot<-
(p_past +
p_future +
plot_layout(guides = "collect",nrow=2, heights=c(1,2))) &
theme(legend.position = "right",strip.background = element_blank(),
strip.text = element_text(face = "bold"))
chill_comp_plotMake a heat map illustrating past and future changes in Safe Winter Chill, relative to a past scenario, for the 13 chill models used here
library(chillR)
##install.packages("remotes")
##install.packages("colorRamps")
install.packages("gganimate")
library(remotes)
remotes::install_github("EduardoFernandezC/dormancyR")
library(dormancyR)
library(ggplot2)
library(reshape2)
library(colorRamps)
library(kableExtra)
library(patchwork)
library(gganimate)
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
hourly_models <- list(Chilling_units = chilling_units,
Low_chill = low_chill_model,
Modified_Utah = modified_utah_model,
North_Carolina = north_carolina_model,
Positive_Utah = positive_utah_model,
Chilling_Hours = Chilling_Hours,
Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model)
daily_models<-list(Rate_of_Chill = rate_of_chill,
Chill_Days = chill_days,
Exponential_Chill = exponential_chill,
Triangula_Chill_Haninnen = triangular_chill_1,
Triangular_Chill_Legave = triangular_chill_2)
metrics<-c(names(daily_models),names(hourly_models))
model_labels=c("Rate of Chill",
"Chill Days",
"Exponential Chill",
"Triangular Chill (Häninnen)",
"Triangular Chill (Legave)",
"Chilling Units",
"Low-Chill Chill Units",
"Modified Utah Chill Units",
"North Carolina Chill Units",
"Positive Utah Chill Units",
"Chilling Hours",
"Utah Chill Units",
"Chill Portions")
kable(data.frame(Metric=model_labels,'Function name'=metrics)) %>%
kable_styling("striped", position = "left",font_size = 10)
MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
Start_JDay<-305
End_JDay<-59
daily_models_past_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
hourly_models_past_scenarios<-tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
#past_scenarios<-daily_models_past_scenarios
past_scenarios<-lapply(
names(daily_models_past_scenarios),
function(x)
cbind(daily_models_past_scenarios[[x]],
hourly_models_past_scenarios[[x]][,names(hourly_models)]))
names(past_scenarios)<-names(daily_models_past_scenarios)
daily_models_observed<-tempResponse_daily(
MApfelstetten_temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_observed<-
daily_models_observed[which(daily_models_observed$Perc_complete>90),]
hourly_models_observed<-tempResponse_daily_list(
MApfelstetten_temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
past_observed<-cbind(
daily_models_observed,
hourly_models_observed[[1]][,names(hourly_models)])
save_temperature_scenarios(past_scenarios,
"data/chill",
"MApfelstetten_multichill_historic")
write.csv(past_observed,
"data/chill/MApfelstetten_multichill_observed.csv",
row.names=FALSE)
##same procedure for future data
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{
Temps<-load_temperature_scenarios(
"data/Weather",
paste0("MApfelstetten_",Time,"_",RCP))
daily_models_future_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_future_scenarios<-lapply(
daily_models_future_scenarios,
function(x) x[which(x$Perc_complete>90),])
hourly_models_future_scenarios<-
tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
future_scenarios<-daily_models_future_scenarios
future_scenarios<-lapply(
names(future_scenarios),
function(x)
cbind(future_scenarios[[x]],
hourly_models_future_scenarios[[x]][,names(hourly_models)]))
names(future_scenarios)<-names(daily_models_future_scenarios)
chill<-future_scenarios
save_temperature_scenarios(
chill,
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
}
##loading and producing chill scenarios
chill_past_scenarios<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_multichill_historic")
chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
chills <-make_climate_scenario(chill_past_scenarios,
caption = "Historic",
historic_data = chill_observed,
time_series = TRUE)
for(RCP in RCPs)
for(Time in Times)
{
chill<-load_temperature_scenarios(
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
if(RCP=="rcp45") RCPcaption <- "RCP4.5"
if(RCP=="rcp85") RCPcaption <- "RCP8.5"
if(Time=="2050") Time_caption <- "2050"
if(Time=="2090") Time_caption <- "2090"
chills <-make_climate_scenario(chill,
caption =c(RCPcaption,Time_caption),
add_to = chills)
}
##heat map for Safe Winter chill
for(i in 1:length(chills))
{ch<-chills[[i]]
if(ch$caption[1]=="Historic")
{GCMs<-rep("none",length(names(ch$data)))
RCPs<-rep("none",length(names(ch$data)))
Years<-as.numeric(ch$labels)
Scenario<-rep("Historic",length(names(ch$data)))} else
{GCMs<-names(ch$data)
RCPs<-rep(ch$caption[1],length(names(ch$data)))
Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
Scenario<-rep("Future",length(names(ch$data)))}
for(nam in names(ch$data))
{for(met in metrics)
{temp_res<-data.frame(Metric=met,
GCM=GCMs[which(nam==names(ch$data))],
RCP=RCPs[which(nam==names(ch$data))],
Year=Years[which(nam==names(ch$data))],
Result=quantile(ch$data[[nam]][,met],0.1),
Scenario=Scenario[which(nam==names(ch$data))])
if(i==1&nam==names(ch$data)[1]&met==metrics[1])
results<-temp_res else
results<-rbind(results,temp_res)
}}}
for(met in metrics)
results[which(results$Metric==met),"SWC"]<-
results[which(results$Metric==met),"Result"]/
results[which(results$Metric==met&results$Year==1980),"Result"]-1
rng = range(results$SWC)
p_future<-ggplot(results[which(!results$GCM=="none"),],
aes(GCM, y=factor(Metric, levels=metrics),
fill = SWC)) +
geom_tile()
p_future
p_future <-
p_future +
facet_grid(RCP ~ Year)
p_future
#make text smaller
p_future <-
p_future +
theme_bw(base_size = 15) +
theme(axis.text = element_text(size=8))
p_future
library(colorRamps)
p_future <-
p_future +
scale_fill_gradientn(colours=matlab.like(15),
labels = scales::percent,
limits=rng)
p_future
p_future <-
p_future +
theme(axis.text.x = element_text(angle = 75, hjust = 1, vjust = 1)) +
labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
scale_y_discrete(labels=model_labels) +
ylab("Chill metric")
p_future
#plot of past scenarios
p_past<-
ggplot(results[which(results$GCM=="none"),],
aes(Year, y=factor(Metric, levels=metrics),
fill = SWC)) +
geom_tile()
p_past<-
p_past +
theme_bw(base_size = 12) +
theme(axis.text = element_text(size=8))
p_past<-
p_past +
scale_fill_gradientn(colours=matlab.like(15),
labels = scales::percent,
limits=rng)
p_past<-
p_past +
scale_x_continuous(position = "top")
p_past<-
p_past +
labs(fill = "Change in\nSafe Winter Chill\nsince 1975") +
scale_y_discrete(labels=model_labels) +
ylab("Chill metric")
chill_comp_plot<-
(p_past +
p_future +
plot_layout(guides = "collect",nrow=2, heights=c(1,2))) &
theme(legend.position = "right",strip.background = element_blank(),
strip.text = element_text(face = "bold"))
chill_comp_plot
Produce an animated line plot of your results (summarizing Safe Winter Chill across all the GCMs).
library(chillR)
##install.packages("remotes")
##install.packages("colorRamps")
install.packages("gganimate")
library(remotes)
remotes::install_github("EduardoFernandezC/dormancyR")
library(dormancyR)
library(ggplot2)
library(reshape2)
library(colorRamps)
library(kableExtra)
library(patchwork)
library(gganimate)
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
hourly_models <- list(Chilling_units = chilling_units,
Low_chill = low_chill_model,
Modified_Utah = modified_utah_model,
North_Carolina = north_carolina_model,
Positive_Utah = positive_utah_model,
Chilling_Hours = Chilling_Hours,
Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model)
daily_models<-list(Rate_of_Chill = rate_of_chill,
Chill_Days = chill_days,
Exponential_Chill = exponential_chill,
Triangula_Chill_Haninnen = triangular_chill_1,
Triangular_Chill_Legave = triangular_chill_2)
metrics<-c(names(daily_models),names(hourly_models))
model_labels=c("Rate of Chill",
"Chill Days",
"Exponential Chill",
"Triangular Chill (Häninnen)",
"Triangular Chill (Legave)",
"Chilling Units",
"Low-Chill Chill Units",
"Modified Utah Chill Units",
"North Carolina Chill Units",
"Positive Utah Chill Units",
"Chilling Hours",
"Utah Chill Units",
"Chill Portions")
kable(data.frame(Metric=model_labels,'Function name'=metrics)) %>%
kable_styling("striped", position = "left",font_size = 10)
MApfelstetten_temps<-read_tab("data/MApfelstetten_temps.csv")
Temps<-load_temperature_scenarios("data/Weather","MApfelstetten_historic")
Start_JDay<-305
End_JDay<-59
daily_models_past_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
hourly_models_past_scenarios<-tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
#past_scenarios<-daily_models_past_scenarios
past_scenarios<-lapply(
names(daily_models_past_scenarios),
function(x)
cbind(daily_models_past_scenarios[[x]],
hourly_models_past_scenarios[[x]][,names(hourly_models)]))
names(past_scenarios)<-names(daily_models_past_scenarios)
daily_models_observed<-tempResponse_daily(
MApfelstetten_temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_observed<-
daily_models_observed[which(daily_models_observed$Perc_complete>90),]
hourly_models_observed<-tempResponse_daily_list(
MApfelstetten_temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
past_observed<-cbind(
daily_models_observed,
hourly_models_observed[[1]][,names(hourly_models)])
save_temperature_scenarios(past_scenarios,
"data/chill",
"MApfelstetten_multichill_historic")
write.csv(past_observed,
"data/chill/MApfelstetten_multichill_observed.csv",
row.names=FALSE)
##same procedure for future data
RCPs<-c("rcp45","rcp85")
Times<-c(2050,2090)
for(RCP in RCPs)
for(Time in Times)
{
Temps<-load_temperature_scenarios(
"data/Weather",
paste0("MApfelstetten_",Time,"_",RCP))
daily_models_future_scenarios<-tempResponse_list_daily(
Temps,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=daily_models)
daily_models_future_scenarios<-lapply(
daily_models_future_scenarios,
function(x) x[which(x$Perc_complete>90),])
hourly_models_future_scenarios<-
tempResponse_daily_list(
Temps,
latitude=48.493,
Start_JDay = Start_JDay,
End_JDay = End_JDay,
models=hourly_models,
misstolerance = 10)
future_scenarios<-daily_models_future_scenarios
future_scenarios<-lapply(
names(future_scenarios),
function(x)
cbind(future_scenarios[[x]],
hourly_models_future_scenarios[[x]][,names(hourly_models)]))
names(future_scenarios)<-names(daily_models_future_scenarios)
chill<-future_scenarios
save_temperature_scenarios(
chill,
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
}
##loading and producing chill scenarios
chill_past_scenarios<-load_temperature_scenarios(
"data/chill",
"MApfelstetten_multichill_historic")
chill_observed<-read_tab("data/chill/MApfelstetten_multichill_observed.csv")
chills <-make_climate_scenario(chill_past_scenarios,
caption = "Historic",
historic_data = chill_observed,
time_series = TRUE)
for(RCP in RCPs)
for(Time in Times)
{
chill<-load_temperature_scenarios(
"data/chill",
paste0("MApfelstetten_multichill_",Time,"_",RCP))
if(RCP=="rcp45") RCPcaption <- "RCP4.5"
if(RCP=="rcp85") RCPcaption <- "RCP8.5"
if(Time=="2050") Time_caption <- "2050"
if(Time=="2090") Time_caption <- "2090"
chills <-make_climate_scenario(chill,
caption =c(RCPcaption,Time_caption),
add_to = chills)
}
##heat map for Safe Winter chill
for(i in 1:length(chills))
{ch<-chills[[i]]
if(ch$caption[1]=="Historic")
{GCMs<-rep("none",length(names(ch$data)))
RCPs<-rep("none",length(names(ch$data)))
Years<-as.numeric(ch$labels)
Scenario<-rep("Historic",length(names(ch$data)))} else
{GCMs<-names(ch$data)
RCPs<-rep(ch$caption[1],length(names(ch$data)))
Years<-rep(as.numeric(ch$caption[2]),length(names(ch$data)))
Scenario<-rep("Future",length(names(ch$data)))}
for(nam in names(ch$data))
{for(met in metrics)
{temp_res<-data.frame(Metric=met,
GCM=GCMs[which(nam==names(ch$data))],
RCP=RCPs[which(nam==names(ch$data))],
Year=Years[which(nam==names(ch$data))],
Result=quantile(ch$data[[nam]][,met],0.1),
Scenario=Scenario[which(nam==names(ch$data))])
if(i==1&nam==names(ch$data)[1]&met==metrics[1])
results<-temp_res else
results<-rbind(results,temp_res)
}}}
for(met in metrics)
results[which(results$Metric==met),"SWC"]<-
results[which(results$Metric==met),"Result"]/
results[which(results$Metric==met&results$Year==1980),"Result"]-1
##animated plot creation
hist_results<-results[which(results$GCM=="none"),]
hist_results$RCP<-"RCP4.5"
hist_results_2<-hist_results
hist_results_2$RCP<-"RCP8.5"
hist_results<-rbind(hist_results,hist_results_2)
future_results<-results[which(!results$GCM=="none"),]
GCM_aggregate<-aggregate(
future_results$SWC,
by=list(future_results$Metric,future_results$RCP,future_results$Year),
FUN=mean)
colnames(GCM_aggregate)<-c("Metric","RCP","Year","SWC")
RCP_Time_series<-rbind(hist_results[,c("Metric","RCP","Year","SWC")],
GCM_aggregate)
chill_change_plot<-
ggplot(data=RCP_Time_series,
aes(x=Year,y=SWC,col=factor(Metric,levels=metrics))) +
geom_line(linewidth=1.3) +
facet_wrap(~RCP,nrow=2) +
theme_bw(base_size=12) +
labs(col = "Change in\nSafe Winter Chill\nsince 1975") +
scale_color_discrete(labels=model_labels) +
scale_y_continuous(labels = scales::percent) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold")) +
ylab("Safe Winter Chill")
chill_change_plot
library(gganimate)
p<-chill_change_plot + transition_reveal(Year)
## the following is needed to bypass a writing permission issue - not sure why but
## this appears to be working
install.packages("gifski")
library(gifski)
animate(p, renderer = gifski_renderer())
anim_save("data/chill_comparison_animation.gif", animation = last_animation())
Exercise:
Provide a brief narrative describing what p-hacking is, and why this is a problematic approach to data analysis.
In regression fitting, a p-value lower than 0.05 means that the null hypothesis can be rejected (thus there is a significant difference). p-hacking is when researchers handling huge data try out numerous statistical analyses to find and fit patterns or trends that are ‘significant’. These published “cherry-picked” correlations can mislead future research, impede scientific progress, and misguide people by attributing results with the wrong parameters.
Provide a sketch of your causal understanding of the relationship between temperature and bloom dates.
Bloom dates<-f(chilling, forcing), which are temperature-dependent. Hence, more heating means less chill and more forcing; this leads to earlier bloom. Specifically, due to the increase in temperatures in autumn-winter, the blooming schedules get impaired, a phenomenon supported by the earlier bloom of Japanese pear investigated by Tominaga et al in 2022. The researchers speculate that the lack of chilling (warmer temperatures) caused metabolic and physiological changes leading to earlier bloom dates or ‘flowering disorder’.
What do we need to know to build a process-based model from this?
From the previous lessons, I learned that process-based models are the most complex models where lots of parameters need to be considered. I also remember that there should be a point at which we draw the line because considering too many factors can multiply errors too; thus we need a moderate number of variables. We need to consider chill and heat accumulation to build the model.
Exercise:
Briefly explain why you shouldn’t take the results of a PLS regression analysis between temperature and phenology at face value. What do you need in addition in order to make sense of such outputs?
We must not over-interpret PLS regression analysis (temperature and phenology) results because this cannot prove anything on its own. It is comparable to finding patterns in the walls of the bathroom and interpreting them as meaningful images. To make sense of things, we must interpret trends side by side with concrete, proven, and scientific basis– in this case, the ecological theories. Ecological theories can be used to set expectations, then PLS analysis can be used to support or challenge the ideas.
Replicate the PLS analysis for the Roter Boskoop dataset that you used in a previous lesson.
library(chillR)
library(kableExtra)
Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
Roter_first<-Roter[,1:2]
Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
Roter_first<-make_JDay(Roter_first)
Roter_first<-Roter_first[,c("Pheno_year","JDay")]
colnames(Roter_first)<-c("Year","pheno")
kable(Roter_first[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
KA_temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
KA_temps<-make_JDay(KA_temps)
kable(KA_temps[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
PLS_results<-PLS_pheno(KA_temps,Roter_first)
kable(PLS_results$PLS_summary[1:5,]) %>%
kable_styling("striped", position = "left",font_size = 8)
dir.create("pictures")
plot_PLS(PLS_results,PLS_results_path = "pictures/PLS_out")
library(ggplot2)
PLS_gg<-PLS_results$PLS_summary
PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-
ISOdate(2001,
PLS_gg$Month[which(PLS_gg$JDay<=0)],
PLS_gg$Day[which(PLS_gg$JDay<=0)])
PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
geom_bar(stat='identity',aes(fill=VIP_importance))
VIP_plot <- VIP_plot +
scale_fill_manual(name="VIP",
labels = c("<0.8", ">0.8"),
values = c("FALSE"="grey", "TRUE"="blue")) +
theme_bw(base_size=12) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank() )
VIP_plot
coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
scale_fill_manual(name="Effect direction",
labels = c("Advancing", "Unimportant","Delaying"),
values = c("-1"="red", "0"="grey","1"="dark green")) +
theme_bw(base_size=12) +
ylab("PLS coefficient") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank() )
coeff_plot
temp_plot<- ggplot(PLS_gg) +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev,ymax=Tmean+Tstdev),
fill="grey") +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==-1),
ymax=Tmean+Tstdev*(VIP_Coeff==-1)),
fill="red") +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==1),
ymax=Tmean+Tstdev*(VIP_Coeff==1)),
fill="dark green") +
geom_line(aes(x=Date,y=Tmean)) +
theme_bw(base_size=12) +
ylab(expression(paste(T[mean]," (°C)")))
temp_plot
library(patchwork)
plot<- (VIP_plot +
coeff_plot +
temp_plot +
plot_layout(ncol=1,
guides = "collect")
) & theme(legend.position = "right",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
ggplot_PLS<-function(PLS_results)
{
library(ggplot2)
PLS_gg<-PLS_results$PLS_summary
PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-
ISOdate(2001,
PLS_gg$Month[which(PLS_gg$JDay<=0)],
PLS_gg$Day[which(PLS_gg$JDay<=0)])
PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
geom_bar(stat='identity',aes(fill=VIP>0.8)) +
scale_fill_manual(name="VIP",
labels = c("<0.8", ">0.8"),
values = c("FALSE"="grey", "TRUE"="blue")) +
theme_bw(base_size=12) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank() )
coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
scale_fill_manual(name="Effect direction",
labels = c("Advancing", "Unimportant","Delaying"),
values = c("-1"="red", "0"="grey","1"="dark green")) +
theme_bw(base_size=12) +
ylab("PLS coefficient") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank() )
temp_plot<- ggplot(PLS_gg) +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev,ymax=Tmean+Tstdev),
fill="grey") +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==-1),
ymax=Tmean+Tstdev*(VIP_Coeff==-1)),
fill="red") +
geom_ribbon(aes(x=Date,ymin=Tmean-Tstdev*(VIP_Coeff==1),
ymax=Tmean+Tstdev*(VIP_Coeff==1)),
fill="dark green") +
geom_line(aes(x=Date,y=Tmean)) +
theme_bw(base_size=12) +
ylab(expression(paste(T[mean]," (°C)")))
library(patchwork)
plot<- (VIP_plot +
coeff_plot +
temp_plot +
plot_layout(ncol=1,
guides = "collect")
) & theme(legend.position = "right",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
}
ggplot_PLS(PLS_results)
Write down your thoughts on why we’re not seeing the temperature response pattern we may have expected. What happened to the chill response?
From the plot, between mid February to April- advancing bloom dates due to high heat accumulation was observed (clear forcing period; almost same trend with the example given for this chapter. Chilling period is not clear, as what was observed in Alexander Lucas pear dataset in Klein-Altendorf. One reason for unclear chilling response is due to cooler climatic conditions–this blurs the relationship between temperature changes and chill accumulation. Klein-Altendorf was one example of a cold location, maybe the Roter-Boskoop data was also captured in a relatively cooler location.
Exercise:
Briefly explain in what climatic settings we can expect PLS regression to detect the chilling phase - and in what settings this probably won’t work.
In settings with relatively warm temperatures (in the case of the given examples, Davis), chilling period can be clearly detected by PLS; in relatively cooler settings (e.g. Klein-Alterdorf, Beijing), chill accumulation and temperature relationship is not monotonic, so results can be fuzzy.
How could we overcome this problem?
Conversion of temperatures to daily rate of chill accumulation. In Chapter 20, this was achieved by using chillR models, particularly with Dynamic Model.
In finding chilling and forcing periods:
-consistent negative coefs
Chilling period-> relationship between bloom dates and daily chill accumulation
Forcing period-> relationship between blooom dates and daily heat accumulation
In delineating, look at results tables, plots, and ecological theory behind it.
Exercise:
Repeat the PLS_chill_force procedure for the ‘Roter
Boskoop’ dataset. Include plots of daily chill and heat
accumulation.
library(chillR)
library(kableExtra)
temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
library(dplyr)
temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
stack_hourly_temps(latitude=48.5)
head(temps_hourly$hourtemps)
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=1,
models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model, GDH = GDH)
)
head(daychill$daily_chill)
dc<-make_daily_chill_plot2(daychill,metrics=c("Chill_Portions"),cumulative=FALSE,
startdate=300,enddate=30,focusyears=c(2008), metriclabels="Chill Portions")
dc<-make_daily_chill_plot2(daychill,metrics=c("Chill_Portions"),cumulative=TRUE,
startdate=300,enddate=30,focusyears=c(2008), metriclabels="Chill Portions")
Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
Roter_first<-Roter[,1:2]
Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
Roter_first<-make_JDay(Roter_first)
Roter_first<-Roter_first[,c("Pheno_year","JDay")]
colnames(Roter_first)<-c("Year","pheno")
##extracting PLS results table for CP
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = "Chill_Portions",
heat_models = "GDH",
runn_means = 11)
plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")
##extracting PLS results table for CU
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = "Utah_Chill_Units",
heat_models = "GDH",
runn_means = 11)
plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")
##extracting PLS results table for CH
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = "Chilling_Hours",
heat_models = "GDH",
runn_means = 11)
plot_PLS(plscf, PLS_results_path= "data/plots/PLS_outputs_11_day")
Above is the PLS analysis plot (using CP model) applied to the Roter-Boskoop dataset. From here, I can roughly say that chilling period started around early December to February, which is quite similar with the example given in class. Forcing period must be around January to bloom date.
Run PLS_chill_force analyses for all three major
chill models. Delineate your best estimates of chilling and forcing
phases for all of them.
library(chillR)
library(kableExtra)
library(ggplot2)
library(dplyr)
temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
stack_hourly_temps(latitude=48.5)
Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
Roter_first<-Roter[,1:2]
Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
Roter_first<-make_JDay(Roter_first)
Roter_first<-Roter_first[,c("Pheno_year","JDay")]
colnames(Roter_first)<-c("Year","pheno")
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=11,
models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model, GDH = GDH)
)
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
heat_models = c("GDH"))
##creation of function for chilling and forcing period shading for CP, CU, and CH
plot_PLS_chill_force<-function(plscf,
chill_metric="Chill_Portions",
heat_metric="GDH",
chill_label="CP",
heat_label="GDH",
chill_phase=c(-48,62),
heat_phase=c(-5,105.5))
{
PLS_gg<-plscf[[chill_metric]][[heat_metric]]$PLS_summary
PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-ISOdate(2001,PLS_gg$Month[which(PLS_gg$JDay<=0)],PLS_gg$Day[which(PLS_gg$JDay<=0)])
PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
chill_start_date<-ISOdate(2001,12,31)+chill_phase[1]*24*3600
chill_end_date<-ISOdate(2001,12,31)+chill_phase[2]*24*3600
heat_start_date<-ISOdate(2001,12,31)+heat_phase[1]*24*3600
heat_end_date<-ISOdate(2001,12,31)+heat_phase[2]*24*3600
temp_plot<- ggplot(PLS_gg) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev ,
ymax=MetricMean + MetricStdev ),
fill="grey") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev * (VIP_Coeff==-1),
ymax=MetricMean + MetricStdev * (VIP_Coeff==-1)),
fill="red") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev * (VIP_Coeff==1),
ymax=MetricMean + MetricStdev * (VIP_Coeff==1)),
fill="dark green") +
geom_line(aes(x=Date,y=MetricMean )) +
facet_wrap(vars(Type), scales = "free_y",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill=paste0("Chill (",chill_label,")"),Heat=paste0("Heat (",heat_label,")")) )) ) +
ggtitle("Daily chill and heat accumulation rates") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
)
VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_bar(stat='identity',aes(fill=VIP>0.8)) +
facet_wrap(vars(Type), scales="free",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill="VIP for chill",Heat="VIP for heat") )) ) +
scale_y_continuous(limits=c(0,max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
ggtitle("Variable Importance in the Projection (VIP) scores") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
) +
scale_fill_manual(name="VIP",
labels = c("<0.8", ">0.8"),
values = c("FALSE"="grey", "TRUE"="blue")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())
coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
facet_wrap(vars(Type), scales="free",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill="MC for chill",Heat="MC for heat") )) ) +
scale_y_continuous(limits=c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
ggtitle("Model coefficients (MC)") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
) +
scale_fill_manual(name="Effect direction",
labels = c("Advancing", "Unimportant","Delaying"),
values = c("-1"="red", "0"="grey","1"="dark green")) +
ylab("PLS coefficient") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())
library(patchwork)
plot<- (VIP_plot +
coeff_plot +
temp_plot +
plot_layout(ncol=1,
guides = "collect")
) & theme(legend.position = "right",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
}
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=11,
models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model, GDH = GDH)
)
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
heat_models = c("GDH"))
##plot (CH)
plot_PLS_chill_force(plscf,
chill_metric="Chilling_Hours",
heat_metric="GDH",
chill_label="CH",
heat_label="GDH",
chill_phase=c(0,0),
heat_phase=c(0,0))
##plot (CU)
plot_PLS_chill_force(plscf,
chill_metric="Utah_Chill_Units",
heat_metric="GDH",
chill_label="CU",
heat_label="GDH",
chill_phase=c(0,0),
heat_phase=c(0,0))
##plot (CP)
plot_PLS_chill_force(plscf,
chill_metric="Chill_Portions",
heat_metric="GDH",
chill_label="CP",
heat_label="GDH",
chill_phase=c(0,0),
heat_phase=c(0,0))For Chilling Hours model, the chilling period is not clear at
all. The model coefficients are unimportant most of the time, so I
cannot delineate the chilling period for this. For the
forcing period, it is most likely around January to April. I consulted
the table of results for the dates, and I arrived at JDay 8
(January 8) to JDay 105.5 (April 15).
For Utah Chilling Units, chilling period is most likely around JDays -15 to 66. Forcing oeriod is most likely from JDays 10 to 105.5.
For Chilling Portions model, chilling period is most likely around JDays -25 to 62, although I must say that the relationship is not that strong here. Forcing period is most likely from JDays 8 to 105.5; this, in contrast with the former, is very easy to spot.
Plot results for all three analyses, including shaded plot areas for the chilling and forcing periods you estimated.
library(chillR)
library(kableExtra)
library(ggplot2)
library(dplyr)
temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
temps_hourly<-stack_hourly_temps(temps,latitude=48.5)
temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") |>
stack_hourly_temps(latitude=48.5)
Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
Roter_first<-Roter[,1:2]
Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
Roter_first<-make_JDay(Roter_first)
Roter_first<-Roter_first[,c("Pheno_year","JDay")]
colnames(Roter_first)<-c("Year","pheno")
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=11,
models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model, GDH = GDH)
)
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
heat_models = c("GDH"))
##creation of function for chilling and forcing period shading for CP, CU, and CH
plot_PLS_chill_force<-function(plscf,
chill_metric="Chill_Portions",
heat_metric="GDH",
chill_label="CP",
heat_label="GDH",
chill_phase=c(-48,62),
heat_phase=c(-5,105.5))
{
PLS_gg<-plscf[[chill_metric]][[heat_metric]]$PLS_summary
PLS_gg[,"Month"]<-trunc(PLS_gg$Date/100)
PLS_gg[,"Day"]<-PLS_gg$Date-PLS_gg$Month*100
PLS_gg[,"Date"]<-ISOdate(2002,PLS_gg$Month,PLS_gg$Day)
PLS_gg[which(PLS_gg$JDay<=0),"Date"]<-ISOdate(2001,PLS_gg$Month[which(PLS_gg$JDay<=0)],PLS_gg$Day[which(PLS_gg$JDay<=0)])
PLS_gg[,"VIP_importance"]<-PLS_gg$VIP>=0.8
PLS_gg[,"VIP_Coeff"]<-factor(sign(PLS_gg$Coef)*PLS_gg$VIP_importance)
chill_start_date<-ISOdate(2001,12,31)+chill_phase[1]*24*3600
chill_end_date<-ISOdate(2001,12,31)+chill_phase[2]*24*3600
heat_start_date<-ISOdate(2001,12,31)+heat_phase[1]*24*3600
heat_end_date<-ISOdate(2001,12,31)+heat_phase[2]*24*3600
temp_plot<- ggplot(PLS_gg) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev ,
ymax=MetricMean + MetricStdev ),
fill="grey") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev * (VIP_Coeff==-1),
ymax=MetricMean + MetricStdev * (VIP_Coeff==-1)),
fill="red") +
geom_ribbon(aes(x=Date,
ymin=MetricMean - MetricStdev * (VIP_Coeff==1),
ymax=MetricMean + MetricStdev * (VIP_Coeff==1)),
fill="dark green") +
geom_line(aes(x=Date,y=MetricMean )) +
facet_wrap(vars(Type), scales = "free_y",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill=paste0("Chill (",chill_label,")"),Heat=paste0("Heat (",heat_label,")")) )) ) +
ggtitle("Daily chill and heat accumulation rates") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
)
VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_bar(stat='identity',aes(fill=VIP>0.8)) +
facet_wrap(vars(Type), scales="free",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill="VIP for chill",Heat="VIP for heat") )) ) +
scale_y_continuous(limits=c(0,max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
ggtitle("Variable Importance in the Projection (VIP) scores") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
) +
scale_fill_manual(name="VIP",
labels = c("<0.8", ">0.8"),
values = c("FALSE"="grey", "TRUE"="blue")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())
coeff_plot<- ggplot(PLS_gg,aes(x=Date,y=Coef)) +
annotate("rect",
xmin = chill_start_date,
xmax = chill_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "blue") +
annotate("rect",
xmin = heat_start_date,
xmax = heat_end_date,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "red") +
annotate("rect",
xmin = ISOdate(2001,12,31) + min(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
xmax = ISOdate(2001,12,31) + max(plscf$pheno$pheno,na.rm=TRUE)*24*3600,
ymin = -Inf,
ymax = Inf,
alpha = .1,fill = "black") +
geom_vline(xintercept = ISOdate(2001,12,31) + median(plscf$pheno$pheno,na.rm=TRUE)*24*3600, linetype = "dashed") +
geom_bar(stat='identity',aes(fill=VIP_Coeff)) +
facet_wrap(vars(Type), scales="free",
strip.position="left",
labeller = labeller(Type = as_labeller(c(Chill="MC for chill",Heat="MC for heat") )) ) +
scale_y_continuous(limits=c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
ggtitle("Model coefficients (MC)") +
theme_bw(base_size=15) +
theme(strip.background = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(size =12),
plot.title = element_text(hjust = 0.5),
axis.title.y=element_blank()
) +
scale_fill_manual(name="Effect direction",
labels = c("Advancing", "Unimportant","Delaying"),
values = c("-1"="red", "0"="grey","1"="dark green")) +
ylab("PLS coefficient") +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank())
library(patchwork)
plot<- (VIP_plot +
coeff_plot +
temp_plot +
plot_layout(ncol=1,
guides = "collect")
) & theme(legend.position = "right",
legend.text = element_text(size=8),
legend.title = element_text(size=10),
axis.title.x=element_blank())
plot
}
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=11,
models = list(Chilling_Hours = Chilling_Hours, Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model, GDH = GDH)
)
plscf<-PLS_chill_force(daily_chill_obj=daychill,
bio_data_frame=Roter_first,
split_month=6,
chill_models = c("Chilling_Hours", "Utah_Chill_Units", "Chill_Portions"),
heat_models = c("GDH"))
##shaded plot (CH)
plot_PLS_chill_force(plscf,
chill_metric="Chilling_Hours",
heat_metric="GDH",
chill_label="CH",
heat_label="GDH",
chill_phase=c(0,0),
heat_phase=c(8,105.5))
##shaded plot (CU)
plot_PLS_chill_force(plscf,
chill_metric="Utah_Chill_Units",
heat_metric="GDH",
chill_label="CU",
heat_label="GDH",
chill_phase=c(-15,66),
heat_phase=c(10,105.5))
##shaded plot (CP)
plot_PLS_chill_force(plscf,
chill_metric="Chill_Portions",
heat_metric="GDH",
chill_label="CP",
heat_label="GDH",
chill_phase=c(-25,62),
heat_phase=c(8,105.5))
For CH model:
For CU model:
For CP model:
Exercise:
Look across all the PLS results presented above. Can you detect a pattern in where chilling and forcing periods could be delineated clearly, and where this attempt failed?
In cold conditions (Beijing), it is very hard to delineate chilling period from the graphs, but forcing period is clear.
In slightly warmer conditions like Klein-Altendorf, it was still difficult to see the chilling period still.
In warmer conditions, delineating chilling period may be easy due to very strong relationships (high VIP and consistently very low model coefficients), they seem to not agree with expected range of dates (extends beyond normal). As conditions get warmer, chilling periods get easier to spot and get longer, while forcing gets harder to delineate.
Based on this, PLS may not be that useful in extreme conditions (very cold and very warm). We saw that the graphs gave confusing trends and sometimes outrageous results compared to expected. Chilling periods are hard to see in colder conditions, and forcing periods are hard to delineate in warmer settings.
Can you think about possible reasons for the success or failure of PLS analysis based on agroclimatic metrics (if so, write down your thoughts)?
In moderate climates (adequate heat and chill during the right time), the chilling and forcing pperiods are easy to delineate using PLS analysis. This might be because they just get enough chill and hit to fulfill in their development stages, so the trends look clear (hits dates we expect based on theory). I think extreme cold and warm settings tend to make things uninterpretable for PLS analysis because there are many possible significant processes, effects ,and relationships arising which are not considered in the models.
Exercise:
Produce chill and heat model sensitivity plots for the location you focused on in previous exercises.
library(chillR)
library(dplyr)
library(lubridate)
library(ggplot2)
library(colorRamps)
library(lubridate) # to get convenient date functionality
library(purrr) # to get the pluck function
##function for chill model sensitivity
Chill_model_sensitivity<-
function(latitude,
temp_models=list(Dynamic_Model=Dynamic_Model,GDH=GDH),
month_range=c(10,11,12,1,2,3),
Tmins=c(-10:20),
Tmaxs=c(-5:30))
{
metrics<-as.list(rep(NA,length(temp_models)))
names(metrics)<-names(temp_models)
first_run<-TRUE
for(mon in month_range)
{weather<-make_all_day_table(
data.frame(Year=c(2001,2001),
Month=c(mon,mon),
Day=c(1,days_in_month(mon)),
Tmin=c(0,0),
Tmax=c(0,0)))
for(tmin in Tmins)
for(tmax in Tmaxs)
if(tmax>=tmin)
{
hourtemps <- weather %>%
mutate(Tmin=tmin, Tmax=tmax) %>%
stack_hourly_temps(latitude=latitude) %>%
pluck(1,"Temp")
metric <- lapply(temp_models,
function(x) do.call(x,list(hourtemps)) %>%
tail(1)/(length(hourtemps)/24))
if(first_run)
{for(tm in 1:length(temp_models))
metrics[[tm]]<-metric[[tm]]
mins<-tmin
maxs<-tmax
month<-mon
first_run<-FALSE} else
{for(tm in 1:length(temp_models))
metrics[[tm]]<-append(metrics[[tm]],metric[[tm]])
mins<-append(mins,tmin)
maxs<-append(maxs,tmax)
month<-append(month,mon)}
}
}
results<-cbind(data.frame(Month=month,Tmin=mins,Tmax=maxs),
as.data.frame(metrics))
results
}
Chill_sensitivity_temps<-function(chill_model_sensitivity_table,
temperatures,
temp_model,
month_range=c(10,11,12,1,2,3),
Tmins=c(-10:20),
Tmaxs=c(-5:30),
legend_label="Chill/day (CP)")
{
library(ggplot2)
library(colorRamps)
cmst<-filter(chill_model_sensitivity_table,Month %in% month_range) %>%
mutate(Month_names = factor(Month,
levels=month_range,
labels=month.name[month_range]))
DM_sensitivity<-
ggplot(cmst,
aes(x=Tmin,y=Tmax,fill=.data[[temp_model]])) +
geom_tile() +
scale_fill_gradientn(colours=alpha(matlab.like(15), alpha = .5),
name=legend_label) +
coord_cartesian(xlim=range(Tmins),
ylim=range(Tmaxs))
temperatures<-temperatures %>%
filter(Month %in% month_range) %>%
filter(Tmax>=Tmin) %>%
mutate(Month_names=factor(Month,
levels=month_range,
labels=month.name[month_range]) )
DM_sensitivity +
geom_point(data=temperatures,
aes(x=Tmin,y=Tmax,fill=NULL,color="Temperature"),
size=0.2) +
facet_wrap(vars(Month_names)) +
scale_color_manual(values = "black",
labels = "Daily temperature \nextremes (°C)",
name="Observed at site" ) +
guides(fill = guide_colorbar(order = 1),
color = guide_legend(order = 2)) +
ylab("Tmax (°C)") +
xlab("Tmin (°C)") +
theme_bw(base_size=12)
}
Model_sensitivities_MApfelstetten<-
Chill_model_sensitivity(latitude=48.5,
temp_models=list(Dynamic_Model=Dynamic_Model,
GDH=GDH),
month_range=c(10:12,1:5))
write.csv(Model_sensitivities_MApfelstetten,
"data/Model_sensitivities_MApfelstetten.csv",row.names = FALSE)
Model_sensitivities_MApfelstetten<-read.csv("data/Model_sensitivities_MApfelstetten.csv")
MApfelstetten_weather<-read_tab("data/MApfelstetten_weather.csv")
Chill_sensitivity_temps(Model_sensitivities_MApfelstetten,
MApfelstetten_weather,
temp_model="Dynamic_Model",
month_range=c(10,11,12,1,2,3),
legend_label="Chill per day \n(Chill Portions)") +
ggtitle("Chill model sensitivity near Münsingen Apfelstetten, Germany")
Chill_sensitivity_temps(Model_sensitivities_MApfelstetten,
MApfelstetten_weather,
temp_model="GDH",
month_range=c(12,1:5),
legend_label="Heat per day \n(GDH)") +
ggtitle("Heat model sensitivity near Münsingen Apfelstetten, Germany")
January and February show the most variation which makes it suitable for PLS analysis, but all the other months seem to have most days within the optimal chilling conditions.
The GDH model sensitivity plot above shows good variation which makes it clear why it was very easy to spot the forcing period in the PLS analysis.
Look at the plots for the agroclimate-based PLS analyses of the ‘Alexander Lucas’ and ‘Roter Boskoop’ datasets. Provide your best estimates of the chilling and forcing phases.
Pasting here the PLS Analysis plot for ‘Alexander Lucas’ using the CP model (from the Tree Phenology using R module Chapter 20.4):
I also looked at the PLS analysis results here, and I think the chilling period is around JDay -31 to JDay 59. Forcing period is most probably from JDay 1 to 105.5 (bloom date).
Pasting here the shaded PLS Analysis plot for ‘Roter Boskoop’ using the CP model from my work for Chapter 20:
Chilling period is most likely around JDays -25 to 62, although I must say that the relationship is not that strong here. Forcing period is most likely from JDays 8 to 105.5; this, in contrast with the former, is very easy to spot.
Exercise:
Reproduce the analysis for the ‘Roter Boskoop’ dataset.
library(chillR)
library(dplyr)
library(purrr)
Roter_first<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")[,1:2] %>%
mutate(Year=substr(First_bloom,1,4),
Month=substr(First_bloom,5,6),
Day=substr(First_bloom,7,8)) %>%
make_JDay() %>%
select(Pheno_year,JDay) %>%
rename(Year=1,pheno=2)
temps<-read_tab("data/TMaxTMin1958-2019_patched.csv")
temps_hourly<-read_tab("data/TMaxTMin1958-2019_patched.csv") %>%
stack_hourly_temps(latitude=50.6)
daychill<-daily_chill(hourtemps=temps_hourly,
running_mean=1,
models = list(Chilling_Hours = Chilling_Hours,
Utah_Chill_Units = Utah_Model,
Chill_Portions = Dynamic_Model,
GDH = GDH)
)
library(ggplot2)
chill_phase<-c(340,62)
heat_phase<-c(357,105.5)
chill<-tempResponse(hourtemps = temps_hourly,
Start_JDay = chill_phase[1],
End_JDay = chill_phase[2],
models = list(Chill_Portions = Dynamic_Model),
misstolerance=10)
heat<-tempResponse(hourtemps = temps_hourly,
Start_JDay = heat_phase[1],
End_JDay = heat_phase[2],
models = list(GDH = GDH))
##plotting
ggplot(data=chill,aes(x=Chill_Portions)) +
geom_histogram() +
ggtitle("Chill accumulation during endodormancy (Chill Portions)") +
xlab("Chill accumulation (Chill Portions)") +
ylab("Frequency between 1958 and 2019") +
theme_bw(base_size = 12)
ggplot(data=heat,aes(x=GDH)) +
geom_histogram() +
ggtitle("Heat accumulation during ecodormancy (GDH)") +
xlab("Heat accumulation (Growing Degree Hours)") +
ylab("Frequency between 1958 and 2019") +
theme_bw(base_size = 12)
chill_requirement <- mean(chill$Chill_Portions) ##[1] 56.16143
chill_req_error <- sd(chill$Chill_Portions) ##[1] 7.126878
median(chill$Chill_Portions) ##[1] 57.6239
heat_requirement <- mean(heat$GDH) ##[1] 3564.06
heat_req_error <- sd(heat$GDH) ##[1] 1428.505
##function creation for the phenological trend
pheno_trend_ggplot<-function(temps,
pheno,
chill_phase,
heat_phase,
phenology_stage="Bloom")
{
library(fields)
library(reshape2)
library(metR)
library(ggplot2)
library(colorRamps)
# first, a sub-function (function defined within a function) to
# compute the temperature means
mean_temp_period<-function(temps,
start_JDay,
end_JDay,
end_season = end_JDay)
{ temps_JDay<-make_JDay(temps)
temps_JDay[,"Season"]<-temps_JDay$Year
if(start_JDay>end_season)
temps_JDay$Season[which(temps_JDay$JDay>=start_JDay)]<-
temps_JDay$Year[which(temps_JDay$JDay>=start_JDay)]+1
if(start_JDay>end_season)
sub_temps<-subset(temps_JDay,JDay<=end_JDay|JDay>=start_JDay)
if(start_JDay<=end_JDay)
sub_temps<-subset(temps_JDay,JDay<=end_JDay&JDay>=start_JDay)
mean_temps<-aggregate(sub_temps[,c("Tmin","Tmax")],
by=list(sub_temps$Season),
FUN=function(x) mean(x, na.rm=TRUE))
mean_temps[,"n_days"]<-aggregate(sub_temps[,"Tmin"],
by=list(sub_temps$Season),
FUN=length)[,2]
mean_temps[,"Tmean"]<-(mean_temps$Tmin+mean_temps$Tmax)/2
mean_temps<-mean_temps[,c(1,4,2,3,5)]
colnames(mean_temps)[1]<-"End_year"
return(mean_temps)
}
mean_temp_chill<-mean_temp_period(temps = temps,
start_JDay = chill_phase[1],
end_JDay = chill_phase[2],
end_season = heat_phase[2])
mean_temp_heat<-mean_temp_period(temps = temps,
start_JDay = heat_phase[1],
end_JDay = heat_phase[2],
end_season = heat_phase[2])
mean_temp_chill<-
mean_temp_chill[which(mean_temp_chill$n_days >=
max(mean_temp_chill$n_days)-1),]
mean_temp_heat<-
mean_temp_heat[which(mean_temp_heat$n_days >=
max(mean_temp_heat$n_days)-1),]
mean_chill<-mean_temp_chill[,c("End_year","Tmean")]
colnames(mean_chill)[2]<-"Tmean_chill"
mean_heat<-mean_temp_heat[,c("End_year","Tmean")]
colnames(mean_heat)[2]<-"Tmean_heat"
phase_Tmeans<-merge(mean_chill,mean_heat, by="End_year")
colnames(pheno)<-c("End_year","pheno")
Tmeans_pheno<-merge(phase_Tmeans,pheno, by="End_year")
# Kriging interpolation
k<-Krig(x=as.matrix(Tmeans_pheno[,c("Tmean_chill","Tmean_heat")]),
Y=Tmeans_pheno$pheno)
pred<-predictSurface(k)
colnames(pred$z)<-pred$y
rownames(pred$z)<-pred$x
melted<-melt(pred$z)
colnames(melted)<-c("Tmean_chill","Tmean_heat","value")
ggplot(melted,aes(x=Tmean_chill,y=Tmean_heat,z=value)) +
geom_contour_fill(bins=60) +
scale_fill_gradientn(colours=alpha(matlab.like(15)),
name=paste(phenology_stage,"date \n(day of the year)")) +
geom_contour(col="black") +
geom_text_contour(stroke = 0.2) +
geom_point(data=Tmeans_pheno,
aes(x=Tmean_chill,y=Tmean_heat,z=NULL),
size=0.7) +
ylab(expression("Forcing phase "*T[mean]*" ("*degree*C*")")) +
xlab(expression("Chilling phase "*T[mean]*" ("*degree*C*")")) +
theme_bw(base_size=15)
}
chill_phase<-c(340,62)
heat_phase<-c(357,105.5)
pheno_trend_ggplot(temps=temps,
pheno=Roter_first,
chill_phase=chill_phase,
heat_phase=heat_phase,
phenology_stage="Bloom")
I will not repeat the creation of the shaded PLS analysis plot I did for ‘Roter Boskoop’ dataset, I just paste it again here:
From the analysis, the following values were obtained: chilling
requirement (chill_requirement) is around 56.2
CP with an estimated error (chill_req_error) of
7.1 CP. The heat requirement
(heat_requirement) is approximately 3564
GDH, with an estimated error (heat_req_error) of
1428 GDH. In temperate regions, chilling requirements
range from 64.7 to 79.4 CP, and heat requirements are
from 3292.0–8383.6 GDH (Xu et al., 2020); which means
the Roter Boskoop’s chilling requirement is lower than normally required
by temperate trees, while its heat requirement is just a bit higher than
the lower bounds. The phenological plot showed very narrow diagonal
shape (strong temperature correlations during eco and endo-dormancy
stages), this one is more narrow compared to the ‘Alexander Lucas’ plot.
Warm forcing phase led to speeding up of bloom (as early as JDay 100),
and cool forcing period led to bloom delay of up to 20 days.
We’ve looked at data from a number of locations so far. How would you expect this surface plot to look like in Beijing? And how should it look in Tunisia?
Beijing has very cold winters, so probably it would just be an extreme, more pronounced version of the thin, diagonal Klein-Altendorf plot. I would be expecting that cool forcing phase would delay bloom (forcing effect is more dominant in cold conditions), and with warm forcing phase comes earlier bloom. Blue would show up in the upper right side and reds will be on the lower left side.
Tunisia is quite warm, and one example already given with relatively warm conditions (e.g. Davis), and it had very broad surface plot. Maybe Tunisia’s would be a bit more extremes since it is hotter. A warm chilling phase with cool forcing would delay bloom dates (red on lower right side). The opposite would lead to earlier bloom (blue at top left).
Exercise:
Describe the temperature response hypothesis outlined in this chapter.
In ecodormancy stage, bud breaks when target heat accumulation is reached. Endodormancy release happens when target chill requirement is reached. The temperature response theory suggests that shorter dormancy arises from exposure to enough chill during endodormancy stage and enough heat during ecodormancy. Summary:
Exercise:
Explain the difference between output validation and process validation.
Output validation only compares (RMSE, Model Efficiency, etc.) the outputs of the models vs observed values; hence neglecting the process (the how’s and why’s). Process validation meshes models with the relevant processes involved; in the context of tree phenology models, the phenological processes should be the reference when doing process validation.
Explain what a validity domain is and why it is important to consider this whenever we want to use our model to forecast something.
Validity domain is the specific set or range of conditions the model can be used for with acceptable degree of correctness in its predictions.
What is validation for purpose?
It is the evaluation of the aim and use cases of the model prior to creating it.
How can we ensure that our model is suitable for the predictions we want to make?
By checking first if the conditions are within the validity domain of the model. If the conditions are not met, then we should not apply the model as is.
PhenoFlex modelExercise:
Parameterize the PhenoFlex model for `Roter Boskoop’
apples.
I encountered an error here with the ‘phenologyFitter’ function. Apparently it checks that the seasons list is the same number with the number of years in the bloom dates dataset. And upon inspection, ‘Roter Boskoop’ only had data til 2018, so I had to change that in the ‘SeasonsList’.
library(chillR)
library(ggplot2)
library(kableExtra)
CKA_weather<-read_tab("data/TMaxTMin1958-2019_patched.csv")
Roter<-read_tab("data/Roter_Boskoop_bloom_1958_2019.csv")
Roter_first<-Roter[,1:2]
Roter_first[,"Year"]<-substr(Roter_first$First_bloom,1,4)
Roter_first[,"Month"]<-substr(Roter_first$First_bloom,5,6)
Roter_first[,"Day"]<-substr(Roter_first$First_bloom,7,8)
Roter_first<-make_JDay(Roter_first)
Roter_first<-Roter_first[,c("Pheno_year","JDay")]
colnames(Roter_first)<-c("Year","pheno")
hourtemps <- stack_hourly_temps(CKA_weather, latitude=50.6)
# here's the order of the parameters (from the helpfile of the
# PhenoFlex_GDHwrapper function)
# yc, zc, s1, Tu, E0, E1, A0, A1, Tf, Tc, Tb, slope
par <- c(40, 190, 0.5, 25, 3372.8, 9900.3, 6319.5, 5.939917e13, 4, 36, 4, 1.60)
upper <- c(41, 200, 1.0, 30, 4000.0, 10000.0, 7000.0, 6.e13, 10, 40, 10, 50.00)
lower <- c(38, 180, 0.1, 0 , 3000.0, 9000.0, 6000.0, 5.e13, 0, 0, 0, 0.05)
SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
Fit_res <- phenologyFitter(par.guess=par,
modelfn = PhenoFlex_GDHwrapper,
bloomJDays=Roter_first$pheno[which(Roter_first$Year>1958)],
SeasonList=SeasonList,
lower=lower,
upper=upper,
control=list(smooth=FALSE, verbose=FALSE, maxit=1000,
nb.stop.improvement=5))
Roter_par<-Fit_res$par
write.csv(Roter_par,"data/PhenoFlex_parameters_Roter_Boskoop.csv")
|
V1 |
V2 |
|---|---|
|
NA |
,“x” |
|
1 |
,40.5138147244476 |
|
2 |
,199.722028620816 |
|
3 |
,0.18019167088005 |
|
4 |
,24.9705987499286 |
|
5 |
,3372.29197129426 |
Produce plots of predicted vs. observed bloom dates and distribution of prediction errors.
library(ggplot2)
library(tidyr)
require(reshape2)
library(chillR)
Roter_par<-read_tab("data/PhenoFlex_parameters_Roter_Boskoop.csv")[,2]
SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
Roter_PhenoFlex_predictions<-Roter_first[which(Roter_first$Year>1958),]
for(y in 1:length(Roter_PhenoFlex_predictions$Year))
Roter_PhenoFlex_predictions$predicted[y]<-PhenoFlex_GDHwrapper(SeasonList[[y]],Roter_par)
Roter_PhenoFlex_predictions$Error<-
Roter_PhenoFlex_predictions$predicted-Roter_PhenoFlex_predictions$pheno
ggplot(Roter_PhenoFlex_predictions,aes(x=pheno,y=predicted)) +
geom_point() +
geom_abline(intercept=0, slope=1) +
theme_bw(base_size = 15) +
xlab("Observed bloom date (Day of the year)") +
ylab("Predicted bloom date (Day of the year)") +
ggtitle("Predicted vs. observed bloom dates")
ggplot(Roter_PhenoFlex_predictions,aes(Error)) +
geom_histogram() +
ggtitle("Distribution of prediction errors")
Compute the model performance metrics RMSEP, mean error and mean absolute error.
library(chillR)
library(ggplot2)
Roter_par<-read_tab("data/PhenoFlex_parameters_Roter_Boskoop.csv")[,2]
SeasonList <- genSeasonList(hourtemps$hourtemps, mrange = c(8, 6), years=c(1959:2018))
Roter_PhenoFlex_predictions<-Roter_first[which(Roter_first$Year>1958),]
for(y in 1:length(Roter_PhenoFlex_predictions$Year))
Roter_PhenoFlex_predictions$predicted[y]<-PhenoFlex_GDHwrapper(SeasonList[[y]],Roter_par)
Roter_PhenoFlex_predictions$Error<-
Roter_PhenoFlex_predictions$predicted-Roter_PhenoFlex_predictions$pheno
RMSEP(Roter_PhenoFlex_predictions$predicted,Roter_PhenoFlex_predictions$pheno)
##[1] 3.869168
RPIQ(Roter_PhenoFlex_predictions$predicted,Roter_PhenoFlex_predictions$pheno)
##[1] 3.682963
mean(Roter_PhenoFlex_predictions$Error)
##[1] -1.077083
mean(abs(Roter_PhenoFlex_predictions$Error))
##[1] 2.859028
RMSEP is 3.87, mean error is -1.08, and mean absolute error is 2.86.
PhenoFlex model - a second lookExercise:
Make chill and heat response plots for the ‘Alexander Lucas’
PhenoFlex model for the location you did the earlier
analyses for.
Since the Phenoflex functions take too much time to load, I did not run the previous code example to get the phenology parameters for ‘Alexander Lucas’ dataset. I merely downloaded from here. I skipped the chill and heat response curves for ‘Alexander Lucas’ dataset because they are already presented in the module.
library(ggplot2)
library(tidyr)
require(reshape2)
library(chillR)
apply_const_temp <- function(temp, A0, A1, E0, E1, Tf, slope, portions=1200, deg_celsius=TRUE)
{
temp_vector <- rep(temp, times=portions)
res <- chillR::DynModel_driver(temp=temp_vector,
A0=A0, A1=A1,
E0=E0, E1=E1,
Tf=Tf,
slope=slope,
deg_celsius=deg_celsius)
return(res$y[length(res$y)])
}
gen_bell <- function(par, temp_values=seq(-5, 20, 0.1)) {
E0 <- par[5]
E1 <- par[6]
A0 <- par[7]
A1 <- par[8]
Tf <- par[9]
slope <- par[12]
y <- c()
for(i in seq_along(temp_values)) {
y[i] <- apply_const_temp(temp=temp_values[i],
A0=A0, A1=A1, E0=E0, E1=E1, Tf=Tf, slope=slope)
}
return(invisible(y))
}
GDH_response<-function(T, par)
{Tb<-par[11]
Tu<-par[4]
Tc<-par[10]
GDH_weight <- rep(0, length(T))
GDH_weight[which(T >= Tb & T <= Tu)] <-
1/2 * (1 + cos(pi + pi * (T[which(T >= Tb & T <= Tu)] - Tb)/(Tu - Tb)))
GDH_weight[which(T > Tu & T <= Tc)] <-
(1 + cos(pi/2 + pi/2 * (T[which(T > Tu & T <= Tc)] -Tu)/(Tc - Tu)))
return(GDH_weight)
}
Alex_par<-read_tab("data/PhenoFlex_parameters_Alexander_Lucas.csv")[,2]
temp_values=seq(-5, 30, 0.1)
temp_response<-data.frame(Temperature=temp_values,
Chill_response=gen_bell(Alex_par, temp_values),
Heat_response=GDH_response(temp_values,Alex_par))
melted_response<-melt(temp_response,id.vars="Temperature")
##for chill and heat efficiency
latitude<-48.5
month_range<-c(10,11,12,1,2,3)
Tmins=c(-20:20)
Tmaxs=c(-15:30)
mins<-NA
maxs<-NA
chill_eff<-NA
heat_eff<-NA
month<-NA
simulation_par<-Alex_par
library(lubridate)
for(mon in month_range)
{
weather<-make_all_day_table(data.frame(Year=c(2001,2001),
Month=c(mon,mon),
Day=c(1,days_in_month(mon)),
Tmin=c(0,0),
Tmax=c(0,0)))
for(tmin in Tmins)
for(tmax in Tmaxs)
if(tmax>=tmin)
{
weather$Tmin<-tmin
weather$Tmax<-tmax
hourtemps<-stack_hourly_temps(weather,
latitude=latitude)$hourtemps$Temp
chill_eff<-c(chill_eff,
tail(PhenoFlex(temp=hourtemps,
times=c(1: length(hourtemps)),
A0=simulation_par[7],
A1=simulation_par[8],
E0=simulation_par[5],
E1=simulation_par[6],
Tf=simulation_par[9],
slope=simulation_par[12],
deg_celsius=TRUE,
basic_output=FALSE)$y,1)/
days_in_month(mon))
heat_eff<-c(heat_eff,
tail(cumsum(GDH_response(hourtemps,
simulation_par)),1)/
days_in_month(mon))
mins<-c(mins,tmin)
maxs<-c(maxs,tmax)
month<-c(month,mon)
}
}
results<-data.frame(Month=month,Tmin=mins,Tmax=maxs,Chill_eff=chill_eff,Heat_eff=heat_eff)
results<-results[!is.na(results$Month),]
write.csv(results,"data/model_sensitivity_PhenoFlex.csv")
Chill_sensitivity_temps<-function(chill_model_sensitivity_table,
temperatures,
temp_model,
month_range=c(10,11,12,1,2,3),
Tmins=c(-20:20),
Tmaxs=c(-15:30),
legend_label="Chill/day (CP)")
{
library(ggplot2)
library(colorRamps)
cmst<-chill_model_sensitivity_table
cmst<-cmst[which(cmst$Month %in% month_range),]
cmst$Month_names<- factor(cmst$Month, levels=month_range,
labels=month.name[month_range])
DM_sensitivity<-ggplot(cmst,aes(x=Tmin,y=Tmax,fill=.data[[temp_model]])) +
geom_tile() +
scale_fill_gradientn(colours=alpha(matlab.like(15), alpha = .5),
name=legend_label) +
xlim(Tmins[1],Tmins[length(Tmins)]) +
ylim(Tmaxs[1],Tmaxs[length(Tmaxs)])
temperatures<-
temperatures[which(temperatures$Month %in% month_range),]
temperatures[which(temperatures$Tmax<temperatures$Tmin),
c("Tmax","Tmin")]<-NA
temperatures$Month_names <- factor(temperatures$Month,
levels=month_range, labels=month.name[month_range])
DM_sensitivity +
geom_point(data=temperatures,
aes(x=Tmin,y=Tmax,fill=NULL,color="Temperature"),
size=0.2) +
facet_wrap(vars(Month_names)) +
scale_color_manual(values = "black",
labels = "Daily temperature \nextremes (°C)",
name="Observed at site" ) +
guides(fill = guide_colorbar(order = 1),
color = guide_legend(order = 2)) +
ylab("Tmax (°C)") +
xlab("Tmin (°C)") +
theme_bw(base_size=15)
}
Model_sensitivities_PhenoFlex<-
read.csv("data/model_sensitivity_PhenoFlex.csv")
MApfelstetten_weather<-read_tab("data/MApfelstetten_weather.csv")
Chill_sensitivity_temps(Model_sensitivities_PhenoFlex,
MApfelstetten_weather,
temp_model="Chill_eff",
month_range=c(10,11,12,1,2,3),
Tmins=c(-20:20),
Tmaxs=c(-15:30),
legend_label="Chill per day \n(arbitrary)") +
ggtitle("PhenoFlex chill efficiency ('Alexander Lucas')")
Chill_sensitivity_temps(Model_sensitivities_PhenoFlex,
MApfelstetten_weather,
temp_model="Heat_eff",
month_range=c(10,11,12,1,2,3),
Tmins=c(-20:20),
Tmaxs=c(-15:30),
legend_label="Heat per day \n(arbitrary)") +
ggtitle("PhenoFlex heat efficiency ('Alexander Lucas')")
PhenoFlex?Exercise:
What was the objective of this work?
To test and recalibrate parameters of Phenoflex under different conditions.
What was the main conclusion?
Phenoflex, like every other model has limitations; it may be that under extreme conditions (very cold and hot temperatures), the model fails to include the necessary parameters taking play in breaking dormancy. More tests are needed.
What experiments could we conduct to test the hypothesis that emerged at the end of the conclusion?
Phenoflex should be tested more with observed data from very cold, very hot climates, and in between (progression), to be able to clearly draw the validity domains.
Exercise:
Download the phenology dataset for the apple cultivar Roter Boskoop from Klein-Altendorf.
|
V1 |
|---|
|
Pheno_year,First_bloom,Full_bloom,Last_bloom |
|
1958,19580509,19580511,19580514 |
|
1959,19590414,19590421,19590504 |
|
1960,19600424,19600505,19600510 |
|
1961,19610409,19610419,19610501 |
|
1962,19620514,19620522,19620530 |
Illustrate the development of the bloom period over the duration of the weather record. Use multiple ways to show this - feel free to be creative.
library(chillR)
library(leaflet)
library(dplyr)
library(reshape2)
library(kableExtra)
library(ggplot2)
library(Kendall)
CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
Roter<-melt(CKA_Roter,
id.vars = "Pheno_year",
value.name="YEARMODA")
Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
Roter<-make_JDay(Roter)
ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
geom_line() +
theme_bw(base_size=15) +
scale_color_discrete(
name="Phenological event",
labels=c("First bloom", "Full bloom", "Last bloom")) +
xlab("Phenological year") +
ylab("Julian date (day of the year)")
ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
geom_line() +
theme_bw(base_size=15) +
scale_color_discrete(name="Phenological event", labels=c("First bloom", "Full bloom", "Last bloom")) +
xlab("Phenological year") +
ylab("Julian date (day of the year)") +
geom_smooth(method = "lm")
ggplot(data=Roter,aes(Pheno_year,JDay,col=variable)) +
geom_smooth() +
theme_bw(base_size=15) +
scale_color_discrete(
name="Phenological event",
labels=c("First bloom", "Full bloom", "Last bloom")) +
xlab("Phenological year") +
ylab("Julian date (day of the year)")
Evaluate the occurrence of frost events at Klein-Altendorf since 1958. Illustrate this in a plot.
library(chillR)
library(leaflet)
library(dplyr)
library(reshape2)
library(kableExtra)
library(ggplot2)
library(Kendall)
CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
frost_df=data.frame(
lower=c(-1000,0),
upper=c(0,1000),
weight=c(1,0))
frost_model<-function(x) step_model(x,frost_df)
hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
frost<-tempResponse(hourly,models=c(frost=frost_model))
frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
by=list(hourly$hourtemps$YEARMODA),
FUN=sum)
Daily_frost<-make_JDay(CKA_weather)
Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
ggplot(data=Daily_frost,aes(Year,JDay,size=Frost_hours)) +
geom_point(col="light blue",alpha = 0.8) +
scale_size(range = c(0, 3),
breaks = c(1, 5, 10, 15, 20),
labels = c("1", "5", "10", "15", "20"),
name="Frost hours") +
theme_bw(base_size=15)
Produce an illustration of the relationship between spring frost events and the bloom period of ‘Roter Boskoop.’
library(chillR)
library(leaflet)
library(dplyr)
library(reshape2)
library(kableExtra)
library(ggplot2)
library(Kendall)
CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
#bloom
Roter<-melt(CKA_Roter,
id.vars = "Pheno_year",
value.name="YEARMODA")
Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
Roter<-make_JDay(Roter)
##frost
frost_df=data.frame(
lower=c(-1000,0),
upper=c(0,1000),
weight=c(1,0))
frost_model<-function(x) step_model(x,frost_df)
hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
frost<-tempResponse(hourly,models=c(frost=frost_model))
frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
by=list(hourly$hourtemps$YEARMODA),
FUN=sum)
Daily_frost<-make_JDay(CKA_weather)
Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
ggplot(data=Roter,
aes(Pheno_year,JDay,col=variable)) +
geom_line() +
theme_bw(base_size=15) +
scale_color_discrete(
name="Phenological event",
labels=c("First bloom", "Full bloom", "Last bloom")) +
xlab("Phenological year") +
ylab("Julian date (day of the year)") +
geom_point(data=Daily_frost,
aes(Year,JDay,size=Frost_hours),
col="light blue",alpha = 0.8) +
scale_size(range = c(0, 3),
breaks = c(1, 5, 10, 15, 20),
labels = c("1", "5", "10", "15", "20"),
name="Frost hours") +
theme_bw(base_size=15)
ggplot(data=Roter,
aes(Pheno_year,JDay,col=variable)) +
geom_line() +
theme_bw(base_size=15) +
scale_color_discrete(
name="Phenological event",
labels=c("First bloom", "Full bloom", "Last bloom")) +
xlab("Phenological year") +
ylab("Julian date (day of the year)") +
geom_point(data=Daily_frost,
aes(Year,JDay,size=Frost_hours),
col="light blue",
alpha = 0.8) +
scale_size(range = c(0, 3),
breaks = c(1, 5, 10, 15, 20),
labels = c("1", "5", "10", "15", "20"),
name="Frost hours") +
theme_bw(base_size=15) +
ylim(c(77,140))
Ribbon_Roter<-dcast(
Roter,Pheno_year ~ variable, value.var = "JDay")
ggplot(data=Ribbon_Roter,aes(Pheno_year)) +
geom_ribbon(aes(ymin = First_bloom, ymax = Last_bloom),
fill = "light gray") +
geom_line(aes(y = Full_bloom)) +
theme_bw(base_size=15) +
xlab("Phenological year") +
ylab("Julian date (day of the year)") +
geom_point(data=Daily_frost,
aes(Year,JDay,size=Frost_hours),
col="light blue",
alpha = 0.8) +
scale_size(range = c(0, 3),
breaks = c(1, 5, 10, 15, 20),
labels = c("1", "5", "10", "15", "20"),
name="Frost hours") +
theme_bw(base_size=15) +
ylim(c(77,140))
# identify frost events that overlap with bloom
lookup_dates<-Ribbon_Roter
row.names(lookup_dates)<-lookup_dates$Pheno_year
Daily_frost[,"First_bloom"]<-
lookup_dates[as.character(Daily_frost$Year),"First_bloom"]
Daily_frost[,"Last_bloom"]<-
lookup_dates[as.character(Daily_frost$Year),"Last_bloom"]
Daily_frost[
which(!is.na(Daily_frost$Frost_hours)),"Bloom_frost"]<-
"Before bloom"
Daily_frost[
which(Daily_frost$JDay>=Daily_frost$First_bloom),"Bloom_frost"]<-
"During bloom"
Daily_frost[
which(Daily_frost$JDay>Daily_frost$Last_bloom),"Bloom_frost"]<-
"After bloom"
Daily_frost[
which(Daily_frost$JDay>180),"Bloom_frost"]<-
"Before bloom"
ggplot(data=Ribbon_Roter,aes(Pheno_year)) +
geom_ribbon(aes(ymin = First_bloom, ymax = Last_bloom),
fill = "light gray") +
geom_line(aes(y = Full_bloom)) +
theme_bw(base_size=15) +
xlab("Phenological year") +
ylab("Julian date (day of the year)") +
geom_point(data=Daily_frost,
aes(Year,JDay,size=Frost_hours,col=Bloom_frost),
alpha = 0.8) +
scale_size(range = c(0, 5),
breaks = c(1, 5, 10, 15, 20),
labels = c("1", "5", "10", "15", "20"),
name="Frost hours") +
scale_color_manual(
breaks=c("Before bloom", "During bloom", "After bloom"),
values=c("light green","red","light blue"),
name="Frost timing") +
theme_bw(base_size=15) +
ylim(c(77,140))
Evaluate how the risk of spring frost for this cultivar has changed over time. Has there been a significant trend?
library(chillR)
library(leaflet)
library(dplyr)
library(reshape2)
library(kableExtra)
library(ggplot2)
library(Kendall)
CKA_Roter<-read.csv("data/Roter_Boskoop_bloom_1958_2019.csv")
CKA_weather<-read.csv("data/TMaxTMin1958-2019_patched.csv")
#bloom
Roter<-melt(CKA_Roter,
id.vars = "Pheno_year",
value.name="YEARMODA")
Roter$Year<-as.numeric(substr(Roter$YEARMODA,1,4))
Roter$Month<-as.numeric(substr(Roter$YEARMODA,5,6))
Roter$Day<-as.numeric(substr(Roter$YEARMODA,7,8))
Roter<-make_JDay(Roter)
##frost
frost_df=data.frame(
lower=c(-1000,0),
upper=c(0,1000),
weight=c(1,0))
frost_model<-function(x) step_model(x,frost_df)
hourly<-stack_hourly_temps(CKA_weather,latitude=50.625)
frost<-tempResponse(hourly,models=c(frost=frost_model))
frost_model_no_summ<-function(x) step_model(x, frost_df, summ=FALSE)
hourly$hourtemps[,"frost"]<-frost_model_no_summ(hourly$hourtemps$Temp)
Daily_frost_hours<-aggregate(hourly$hourtemps$frost,
by=list(hourly$hourtemps$YEARMODA),
FUN=sum)
Daily_frost<-make_JDay(CKA_weather)
Daily_frost[,"Frost_hours"]<-Daily_frost_hours$x
Daily_frost$Frost_hours[which(Daily_frost$Frost_hours==0)]<-NA
# identify frost events that overlap with bloom
lookup_dates<-Ribbon_Roter
row.names(lookup_dates)<-lookup_dates$Pheno_year
Daily_frost[,"First_bloom"]<-
lookup_dates[as.character(Daily_frost$Year),"First_bloom"]
Daily_frost[,"Last_bloom"]<-
lookup_dates[as.character(Daily_frost$Year),"Last_bloom"]
Daily_frost[
which(!is.na(Daily_frost$Frost_hours)),"Bloom_frost"]<-
"Before bloom"
Daily_frost[
which(Daily_frost$JDay>=Daily_frost$First_bloom),"Bloom_frost"]<-
"During bloom"
Daily_frost[
which(Daily_frost$JDay>Daily_frost$Last_bloom),"Bloom_frost"]<-
"After bloom"
Daily_frost[
which(Daily_frost$JDay>180),"Bloom_frost"]<-
"Before bloom"
##checking for bloom trends
Bloom_frost_trend<-aggregate(
Daily_frost$Frost_hours,
by=list(Daily_frost$Year,Daily_frost$Bloom_frost),
FUN=function(x) sum(x,na.rm=TRUE))
colnames(Bloom_frost_trend)<-c("Year","Frost_timing","Frost_hours")
DuringBloom<-
Bloom_frost_trend[which(Bloom_frost_trend$Frost_timing=="During bloom"),]
ggplot(data=DuringBloom,aes(Year,Frost_hours)) +
geom_col() +
ylab("Frost hours")
Kendall(x=DuringBloom$Year,y=DuringBloom$Frost_hours)
##tau = 0.0333, 2-sided pvalue =0.74038
lm(DuringBloom$Frost_hours~DuringBloom$Year)
##Call:
#lm(formula = DuringBloom$Frost_hours ~ DuringBloom$Year)
#Coefficients:
# (Intercept) DuringBloom$Year
# -91.48652 0.04786
Trend above is hard to interpret, I will not attempt to. And after checking via Kendall test, trend is insignificant (p-value >0.05) anyway.
I am glad I enlisted in this course. I learned from bottom and now I
can say I am still at the bottom but a bit with clue what is going on
(in terms of R, climate scenarios, and tree phenology). The errors I
encountered forced me to dissect function per function, objects, lists,
files, and the code-- it was dreadful but necessary.
Upon using ‘chillR’, I observed some things and I list them down below,
hoping the developers will take a look at them for future
improvements:
loading bars with %progress and time estimates- I am not sure if this is entirely possible, but I am suggesting ‘time estimates’ anyway--would be very helpful in deciding whether one should take a nap while waiting or not.
update lecture portion in module (html) concurrently with the little edits done in the exercises portion-since I am one of those not-so-conventional “studious” students who do exercises in pace with class (which I deeply regretted in the end), I had to rely on notes and modules. In some exercises, some datasets were mentioned, and there were phrases like “that you used in previous chapter”, but they were not even mentioned in the previous chapter module. Bit confusing but figured out stuff by revisiting notes and Eike’s github.
in Chapter 14 there is a persistent error with the ‘getClimateWizardData’ function when dealing with 2080 and 2085 specifically, for RCP8.5 for my location (48.4931° N, 9.3998° E) but it works with the other years (2050 and 2090, for example). It shows the error below.
I tried comparing the outputs and all, but my guess is that there is a special range of temperature considered for temperature scenarios that might exclude some datasets/fail to recognize them. Had to improvise and use trial and error for the year. I learned a lot here though.
#'Warning: scenario_type doesn't look right - is this really a relative scenario?Warning: updating scenario_type to 'absolute'Error: temperature_scenario is NULL
patch_daily_temps seemed to only work perfectly fine with
‘gsod’ datafiles. So many hiccups because I could not also clean the
‘gsod’ raw files with handle_gsod for some reason. It
always returns an empty file even when I tested the Bonn weather codes
and samples from class. That was the reason I had to use ‘dwd’ datafiles
and move the location to Swabian Orchard paradise. After so many
attempts, I was able to figure out that the ‘old and deprecated’
function patch_daily_temperatures would work in my
case.With every word, I felt the passion that the module creators had while crafting the modules. It was very well-made to the point that someone clueless like I exited the course with huge progress in terms of R knowledge. I want more exposure to these concepts so I will be seeing the team again next semester for the Decisions Analysis and Forecasting subject. I will strive to do better for sure. Bis bald!
Jaffrés, Jasmine. (2019). GHCN-Daily: a treasure trove of climate data awaiting discovery. Computers & Geosciences. 122. 35-44. 10.1016/j.cageo.2018.07.003.
Glaser, Rüdiger & Himmelsbach, Iso & Boesmeier, Annette. (2017). Climate of migration? How climate triggered migration from southwest Germany to North America during the 19th century. Climate of the Past. 13. 1573-1592. 10.5194/cp-13-1573-2017.
Kaufmann, H. and Blanke, M.M. (2019). Chilling requirements of
Mediterranean fruit crops in a changing climate. Acta Hortic. 1242,
275-280. DOI: 10.17660/ActaHortic.2019.1242.38
https://doi.org/10.17660/ActaHortic.2019.1242.38
Luedeling, Eike, Evan H Girvetz, Mikhail A Semenov, and Patrick H Brown. 2011. "Climate Change Affects Winter Chill for Temperate Fruit and Nut Trees." PloS One 6 (5): e20155. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0020155.
Milech, Chaiane & Scariotto, Silvia & Dini, Maximiliano & Herter, Flavio & Raseira, Maria. (2018). MODELS TO ESTIMATE CHILLING ACCUMULATION UNDER SUBTROPICAL CLIMATIC CONDITIONS IN BRAZIL. Revista Brasileira de Climatologia. 23. 10.5380/abclima.v23i0.53086. Tominaga, A., Ito, A., Sugiura, T., & Yamane, H. (2022). How Is Global Warming Affecting Fruit Tree Blooming? “Flowering (Dormancy) Disorder” in Japanese Pear (Pyrus pyrifolia) as a Case Study. Frontiers in Plant Science, 12. https://doi.org/10.3389/fpls.2021.787638
The Norns (2018) Norse Mythology for Smart People. Available at: https://norse-mythology.org/gods-and-creatures/others/the-norns/#:~:text=In%20Norse%20mythology%2C%20the%20Norns,any%20and%20all%20other%20beings. (Accessed: March 15, 2023).
Tominaga, A. et al. (2022) "How is global warming affecting Fruit Tree Blooming? 'Flowering (dormancy) disorder' in Japanese pear (Pyrus pyrifolia) as a case study," Frontiers in Plant Science, 12. Available at: https://doi.org/10.3389/fpls.2021.787638.
Schwäbisches Streuobstparadies (no date) Streuobstparadies. Available at: https://www.streuobstparadies.de/ (Accessed: March 15, 2023).
Xu, Yunjia & Dai, Junhu & Ge, Quansheng & Wang, Huanjiong & Tao, Zexing. (2020). Comparison of chilling and heat requirements for leaf unfolding in deciduous woody species in temperate and subtropical China. International Journal of Biometeorology. 10.1007/s00484-020-02007-7.